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Nonorthogonal  wavelet  transforms  play  an  important  role  in  signal  processing 
by  offering  much  greater  flexibility  in  the  choice  of  analyzing  wavelet  than  their 
orthogonal  counterparts;  however,  very  little  literature  has  been  devoted  to  their 
discrete  implementation.  In  this  work  we  introduce  a  new  nonorthogonal  wavelet, 
the  sine-Gabor  function,  and  propose  an  invertible  discrete  implementation  of  the 
wavelet  transform  using  this  new  analyzing  wavelet. 

We  show  that  the  sine-Gabor  function  not  only  satisfies  the  wavelet  admissibil- 
ity condition  but  achieves  nearly  optimum  time-frequency  localization  and  generates 
a  frame  of  wavelets.  Furthermore,  we  show  that  there  is  a  trade-off  between  the  tight- 
ness of  the  frame  and  the  time-frequency  localization  properties  of  the  sine-Gabor 
wavelet.  We  demonstrate  how  this  trade-off  can  be  overcome  by  introducing  voices. 
We  also  show  that  under  certain  conditions  the  sine-Gabor  wavelet  is  nearly  equal 
to  the  first  derivative  of  a  Gaussian  and  generates  a  "snug"  frame.  This  implies  that 
the  first  derivative  of  a  Gaussian  (Canny's  approximation  to  the  optimal  step  edge 
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detector)  can  be  used  to  completely  characterize  and  reconstruct  a  signal  from  its 
discrete  wavelet  coefficients.  We  address  significant  technical  issues  regarding  the  dis- 
cretization of  the  wavelet  transform  and  its  computation  when  using  nonorthogonal 
wavelets.  We  show  that  there  exists  a  class  of  wavelets  (the  sine-Gabor  wavelet  being 
one  of  them)  for  which  the  discrete  wavelet  transform  is  not  invertible.  We  describe 
how  to  evade  this  problem  by  modifying  the  strategy  for  selecting  the  discrete  filters 
implementing  the  transform.  We  show  how  in  the  case  of  the  sine-Gabor  wavelet, 
desirable  properties  such  as  constant  phase  and  nearly  optimum  time-frequency  lo- 
calization are  preserved  by  the  alternative  filters.  We  conclude  by  extending  the 
transform  to  two  dimensions  and  demonstrating  its  use  in  image  processing  with 
both  broad-band  and  narrow-band  sine-Gabor  wavelets  and  as  a  multivoice  multires- 
olution  edge  detector. 
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CHAPTER  1 
INTRODUCTION 

The  analysis  of  nonstationary  signals  often  involves  a  compromise  between 
how  well  transitions  can  be  located,  and  how  finely  long-term  behavior  can  be  iden- 
tified [Vet92].  Wavelets  are  rapidly  finding  application  as  a  tool  for  the  analysis  of 
such  signals  [Com90.  Dau92].  However,  research  has  been  primarily  devoted  to  or- 
thogonal and  biorthogonal  wavelet  transforms  because  of  their  efficient  and  exact  dis- 
crete implementations  [Dau88,  Mal89a,  Mal89b,  Mal89c,  Dau90,  Vet92].  Nonorthog- 
onal  wavelet  transforms  play  an  important  role  in  signal  processing  by  offering  much 
greater  flexibility  in  the  choice  of  analyzing  wavelet  than  their  orthogonal  counter- 
parts. In  particular,  they  allow  adjustment  of  the  relative  resolution  in  time  and 
scale  [She96].  However,  with  some  notable  exceptions  [She92.  She93.  She96].  very 
little  literature  has  been  devoted  to  their  discrete  implementation.  In  this  disserta- 
tion we  introduce  a  new  nonorthogonal  wavelet,  the  sine-Gabor  function,  and  extend 
the  results  of  [She92.  She93,  She96]  to  obtain  an  invertible  discrete  implementation 
of  the  wavelet  transform  when  using  this  new  analyzing  wavelet.  Chapter  2  reviews 
the  theory  of  the  continuous  wavelet  transform  and  wavelet  frames,  listing  a  num- 
ber of  technical  results  that  are  used  later  in  our  study.  Chapter  3  reviews  some  of 
the  connections  between  discrete  filters  and  wavelets  with  an  emphasis  on  the  algo- 
rithms proposed  in  [She92,  She93,  She96]  for  the  computation  of  the  discrete  wavelet 
transform  and  its  inverse  for  nonorthogonal  wavelets.  In  Chapter  4  we  introduce 
the  sine-Gabor  wavelet.  We  show  that  this  function  not  only  satisfies  the  wavelet 
admissibility  condition  but  achieves  nearly  optimum  time-frequency  localization  and 
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generates  a  frame  of  wavelets.  Furthermore,  we  show  that  there  is  a  trade-off  be- 
tween the  tightness  of  the  frame  and  the  time-frequency  localization  properties  of  the 
sine-Gabor  wavelet.  We  demonstrate  how  this  trade-off  can  be  overcome  by  intro- 
ducing voices.  We  also  show  that  under  certain  conditions  the  sine-Gabor  wavelet 
is  nearly  equal  to  the  first  derivative  of  a  Gaussian  and  generates  a  "snug"  frame. 
This  implies  that  the  first  derivative  of  a  Gaussian  (Canny's  approximation  to  the 
optimal  step  edge  detector  [Can86])  can  be  used  to  completely  characterize  and  re- 
construct a  signal  from  its  discrete  wavelet  coefficients.  In  Chapter  5  we  address  the 
implementation  issues  regarding  the  computation  of  the  discrete  wavelet  transform. 
We  show  that  there  exists  a  class  of  wavelets  (the  sine-Gabor  wavelet  being  one  of 
them)  for  which  the  discrete  wavelet  transform  reviewed  in  Chapter  3  is  not  invert- 
ible.  We  describe  how  to  evade  this  problem  by  modifying  the  strategy  for  selecting 
the  discrete  filters  implementing  the  transform.  We  show  how  in  the  case  of  the 
sine-Gabor  wavelet,  desirable  properties  such  as  constant  phase  and  nearly  optimum 
time- frequency  localization  are  preserved  by  the  alternative  filters.  In  Chapter  6  we 
investigate  the  application  of  the  sine-Gabor  wavelet  transform  to  image  processing. 
We  describe  an  extension  of  the  transform  to  two  dimensions  and  demonstrate  the 
use  of  the  transform  with  both  broad-band  and  narrow-band  sine-Gabor  wavelets 
having  nearly  optimum  time-frequency  localization.  We  also  revisit  the  use  of  the 
sine-Gabor  wavelet  as  an  approximation  to  the  first  derivative  of  a  Gaussian  and 
demonstrate  its  use  as  a  multivoice  multiresolution  edge  detector.  Finally,  in  Chap- 
ter 7  we  summarize  our  work  and  describe  possible  extensions  and  open  problems 
that  are  worthy  of  investigation. 


CHAPTER  2 
CONTINUOUS  WAVELET  TRANSFORM 


In  this  chapter  we  review  the  continuous  wavelet  transform  as  well  as  wavelet 
frames.  To  highlight  and  contrast  the  advantages  of  the  wavelet  transform  we  also 
overview  the  short-time  Fourier  transform. 


The  goal  of  signal  analysis  is  to  extract  relevant  information  from  a  signal 
by  transforming  it  into  a  representation  where  the  properties  of  the  signal  are  more 
evident.  For  the  analysis  of  stationary  signals,  that  is,  signals  whose  properties  do 
not  evolve  with  time,  one  example  of  such  representation  is  the  Fourier  transform 


The  Fourier  transform  can  be  viewed  as  the  inner  product  of  the  signal  s(t)  and  the 
sinusoidal  wave  eJwt.  The  analysis  coefficient  s(u>)  =  {Fs){uj)  measures  the  "strength" 
of  the  sinusoidal  wave  of  frequency  u  in  the  signal  s(t). 

Fourier  analysis  works  well  for  signals  composed  of  a  few  stationary  compo- 
nents. However,  any  abrupt  change  in  time  in  a  nonstationary  signal  s(t)  is  spread 
out  over  the  whole  frequency  axis  in  s(u).  In  order  to  adapt  the  Fourier  transform  to 
nonstationary  signals,  Gabor  [Gab46]  introduced  a  new  transform  by  using  a  window 
function  w(t)  in  the  Fourier  integral 


2.1    The  Short-Time  Fourier  Transform 


3 
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A  function  w(t)  qualifies  as  a  window  function  if  it  is  possible  to  identify  its  center 
and  standard  deviation  (or  root  mean  square  (RMS)  duration)  defined  by 


i  r°° 

m(w)  =  — —  /  t\w{t)\2dt 

\H\-2  J-oo 


and 

a(w)=*y<"(t-m(w))2\w(t)\2dt}     ,  (2.1) 

respectively  [Chu92].  In  the  original  Gabor  transform,  the  window  function  was  a 
Gaussian  [Gab46],  however,  the  transform  is  valid  for  any  type  of  window  function 
and  is  referred  to  as  the  short-time  or  windowed  Fourier  transform. 

The  short-time  Fourier  transform  can  be  viewed  as  the  inner  product  of  the 
signal  s(t)  with  the  family  of  functions 


w 


that  is, 

(gws)(r,u)  =  (s(t),  wUt)).  (2-2) 

It  is  easy  to  show  that  wT^(t)  is  a  window  function  with  center  and  standard  deviation 
given  by 

m(uviW)  =  m(w)  +  t  (2.3) 

and 

c(»'r,u..)  =  <r{w),  (2-4) 

respectively.  From  Equations  2.2,  2.3  and  2.4  we  have  that  in  the  time  domain 
the  analysis  coefficient  {Qws)(t,uj)  essentially  depends  on  the  values  of  s(t)  for  t  G 
[m(w)  +  t  -  a(w),  m(w)  +  t  +  <r(w)] .  This  is  called  time  localization.  It  follows  from 
the  above  observation  that  two  pulses  in  time  can  be  discriminated  only  if  they  are 
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more  than  2a(w)  apart.  This  is  referred  to  as  the  time  resolution  of  the  short-time 
Fourier  transform. 

By  applying  Parseval's  theorem  to  Equation  2.2  we  obtain  an  alternative  for- 
mula for  the  short-time  Fourier  transform 

(Gws)(t,u)  =  —  (s(7),  wTAlf))  i  (2-5) 

where 

wrA-y)  =  Mi  -  w)e-i(7-w)T,  (r,w)  e  m2. 

Suppose  10(7)  =  (Jit;) (7)  is  also  a  window  function  with  center,  m(w),  and  standard 
deviation,  a(w).  Then  it  is  easy  to  show  that  wT^{l)  is  a  window  function  with 
center  and  standard  deviation  given  by 

m{'WT,u)  —  m(w)  +  u!  (2.6) 

and 

respectively.  From  Equations  2.5,  2.6  and  2.7  we  have  that  in  the  frequency  domain 
the  analysis  coefficient  {Gws)(t,u)  essentially  depends  on  the  values  of  5(7)  for  7  G 
[m(w)  +u  -  cr(w),  m(w)  +  u  +  a{w)\ .  This  is  called  frequency  localization.  It  follows 
from  the  above  observation  that  two  pure  sinusoids  can  be  discriminated  only  if  their 
frequencies  are  more  than  2<r(w)  apart.  This  is  referred  to  as  the  frequency  resolution 
of  the  short-time  Fourier  transform. 

From  Equations  2.2  through  2.7  we  have  that  (Gws)(t,uj)  yields  a  time- 
frequency  representation  of  s(t).  The  analysis  coefficient  (Gws)(t,u>)  depends  on 
the  time-frequency  window 

[m(it')  +  t  —  a(iv),  m(w)  +  r  +  a(w)]  x  [ra(it))  +  u  —  cr(i£'),  m(w)  +  u  +  cr(io)] . 
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Increasingly  accurate  localization  in  time  and  frequency  is  not  possible  because  the 
area  of  the  time-frequency  window  is  lower  bounded  by  2,  that  is 

Aa{w)a(w)  >  2.  (2.8) 

This  is  referred  to  as  the  uncertainty  principle  or  Heisenberg  inequality.  Equality  is 
satisfied  if  and  only  if  the  window  function  is  a  Gaussian  as  in  the  Gabor  transform. 

The  short-time  Fourier  transform  is  an  isometry  (to  a  proportionality  coeffi- 
cient) from  L2(R)  into  L2(K2)  [Mal89c],  that  is 

1  r+oo  r+oc 

^11  "'ll'i  J -oo    J -oo 

The  function  s(t)  is  reconstructed  from  {Gws)(t,u)  with  the  formula 

i  n+oo  r+oo 

s(t)  =  —±—  /       /      (gws)(r,u)w(t  -  r)e^  drcU. 

The  short-time  Fourier  transform  is  a  redundant  representation.  Instead  of 
computing  {Gws){t,lo)  for  all  values  (r,w)  in  E2,  it  is  possible  to  uniformly  sample 
both  r  and  u>  such  that  the  representation  is  complete  and  stable  [Mal89c].  Let  r0 
and  u0  be  the  sampling  intervals  in  the  time  and  frequency  domains  respectively. 
Then,  the  discrete  short-time  Fourier  transform  is  defined  by  [Mal89c] 

{Gws)m,n  =  (Gws)(™>T0,nu)0),  (m,n)  6  Z2. 

To  reconstruct  any  function  s{t)  G  L2(R)  from  the  set  of  samples  {Gws)m,n,  (m,n)  € 
Z2,  the  operator 

L2(R)  %  £2(Z2) 

must  be  invertible  on  its  range  and  have  a  bounded  inverse  [Mal89c].  In  order  to 
invert  GUM  Daubechies  [Dau90]  has  shown  that  r0  and  uj0  must  verify 

luoTq  <  2ir. 
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A  drawback  of  the  short-time  Fourier  transform  is  that  once  a  window  function 
has  been  chosen  the  time-frequency  resolution  of  the  transform  is  fixed  over  the  entire 
time-frequency  plane.  It  is  therefore  impossible  with  the  short-time  Fourier  transform 
to  analyze  at  the  same  time  transient  signal  components  with  good  time  resolution 
and  quasi-stationary  signal  components  with  good  frequency  resolution  [Rio91]. 

2.2    The  Wavelet  Transform 

The  wavelet  transform  overcomes  the  resolution  limitation  of  the  short-time 
Fourier  transform  by  letting  the  time  and  frequency  resolution  vary  in  the  time- 
frequency  plane.  The  wavelet  transform  is  defined  by  [Chu92] 

(VlV0(a,&)  =  M~*  f°° s(t)eL2(R),a.beR.  (a  +  0).  (2.9) 

In  order  to  reconstruct  s(t)  from  its  wavelet  transform,  the  Fourier  transform  of  tp(t) 
must  satisfy  [Chu92] 

C,.r  a*.  (2.10) 

A  function  ip(t)  is  said  to  be  a  basic  wavelet  if  it  is  possible  to  reconstruct  s{t)  from 
its  wavelet  transform. 

For  any  tp(t)  satisfying  Equation  2.10  we  have  that  the  wavelet  transform  is 
an  isometry  (to  a  proportionality  coefficient)  from  L2(R)  into  I2(R2)  [Chu92],  that 
is 

Nil  =  7^/      /      \ms)(a,b)\2  ^db.  (2.11) 

The  function  s(t)  is  reconstructed  from  (Wts)(a,b),  nj  6  I.  {a  /  0),  with  the 
formula  [Chu92] 

1    p+oo  ^+00  r  /t  -b\]  da 


s 


The  wavelet  transform  can  be  viewed  as  the  inner  product  of  the  signal  s(t) 
with  the  family  of  functions 

Mt)  =  \a\~h  (^)  ,  M€  B,  (a  *  0),  (2.12) 

that  is, 

(W,«)(a,6)  =  <«(*).&*(*)>•  (2-13) 

Suppose  is  a  window  function  with  center  m(ip)  and  standard  deviation 
a(tp).  Then,  it  is  easy  to  show  that  i'a-b{t)  is  a  window  function  with  center  and 
standard  deviation  given  by 

m(i>a;b)  =  am(y)  +  b  (2.14) 

and 

a{ipa-,b)  =  |a|a(V),  (2-15) 

respectively.  From  Equations  2.13,  2.14  and  2.15  we  have  that  in  the  time  domain 
the  analysis  coefficient  (W1j,s)(a.b)  essentially  depends  on  the  values  of  s(t)  for  t  G 
[am(v)  +  b-  |a|a(V'),am(^)  +  b  +  \a\a(if>)] .  It  can  be  verified  [Chu92]  that  if  a  basic 
wavelet,  ip(t),  is  also  a  window  function  then  t/)(t)  is  necessarily  in  I'(E)  so  that  its 
Fourier  transform,  i>(uj)  is  a  continuous  function.  It  follows  from  Equation  2.10  that 
ip(u})  must  vanish  at  the  origin,  that  is, 

/oo 
ij>(t)  dt  =  0.  (2.16) 

By  applying  Parseval's  theorem  to  Equation  2.13  we  obtain  an  alternative 
formula  for  the  wavelet  transform 

(WV)s)M)  =  ^(%),4;6(7)),  (2-17) 
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where 

A*(7)  =  |a|M(a7)e"J76,  a,b  €  R  (a  #  0). 

Suppose  ^(7)  is  also  a  window  function  with  center,  m(ip),  and  standard  deviation, 
cr(^).  Then  it  is  easy  to  show  that  tj^bil)  is  a  window  function  with  center  and 
standard  deviation  given  by 

m(^)  =  (2.18) 

and 

respectively.  From  Equations  2.17,  2.18  and  2.19  we  have  that  in  the  frequency 
domain  the  analysis  coefficient  (W^s)(a,ft)  essentially  depends  on  the  values  of  s{f) 

tor  7  e  |oi  .   a       |0j  j . 

From  Equations  2.13  through  2.19  we  have  that  (W*a)(a,&)  yields  a  time- 
frequency  representation  of  The  analysis  coefficient  (W,,s)(a,&)  depends  on  the 
time-frequency  window 

m(^)     cr(?/>)  m(^)  cr(^) 
[arn(V')  +  6  -  |a|<r(^),  am(^)  +  b+  |a|a(^)]  x    —  rp  —  j^p 

In  contrast  to  the  short-time  Fourier  transform,  both  the  time  and  frequency  res- 
olution in  the  wavelet  transform  vary  in  the  time-frequency  plane.  The  following 
properties  derive  from  the  time-frequency  localization  characteristics  of  the  wavelet 
transform: 

1.  The  ratio  of  the  standard  deviation,  a(ip)/\a\.  of  the  frequency  window  and 
its  center  frequency,  m{xp)/a,  is  given  by  ±^|j,  which  is  independent  of  the 
location  of  the  center  frequency.  This  is  called  constant- Q  frequency  analysis. 


LO 


2.  The  time-frequency  window  narrows  in  time  and  widens  in  frequency  for  large 
center  frequency  m(^)/o  (small  a  >  0),  and  widens  in  time  and  narrows  in 
frequency  for  small  center  frequency  m(ip)/a  (large  a  >  0). 

Recall  that  the  wavelet  transform  in  Equation  2.9  can  be  viewed  as  the  inner 
product  of  s(t)  with  the  family  of  functions  ipa-b(t)  defined  in  Equation  2.12.  Notice 
that  this  family  of  functions  is  obtained  from  a  single  basic  wavelet  by  translating  it  by 
b  and  dilating  it  by  a.  The  basic  wavelet  is  usually  referred  to  as  the  mother  wavelet. 
The  constant  |a|~3  in  Equation  2.12  is  used  for  energy  normalization  purposes,  that 
is  ||^a;6(0ll2  =  11^(0 111-  In  general  the  mother  wavelet  ij>(t)  is  normalized  to  have  its 
energy  equal  to  one,  that  is  \\ip{t)\\l  =  1.  We  will  assume  this  restriction  from  now 
on. 

In  practice  the  parameter  a  is  restricted  to  a  >  0,  so  that  the  center  frequency 
of  Vw>(7)  in  Equation  2.18  is  restricted  to  positive  frequencies.  In  this  context  the 
parameter  a  is  usually  referred  to  as  the  scale  parameter.  In  order  to  reconstruct  s{t) 
from  its  wavelet  transform  restricted  to  a  >  0.  the  Fourier  transform  of  ip{t)  must 
satisfy  a  stricter  admissibility  condition  than  Equation  2.10  given  by  [Chu92] 

r  mt =  r  m-M  „ = \c„  <  x.  (2.20 

Jo  7  Jo  7  2 

For  any  ip{t)  satisfying  Equation  2.20  we  have  the  following  reconstruction  formula 

[Chu92] 

The  wavelet  transform  is  a  redundant  representation.  Instead  of  computing 
(W1f,s)(a,  b)  for  6  e  R  and  a  G  R+,  it  is  possible  to  choose  an  exponential  sampling  of 
the  scale  parameter,  a  =  a£,  j  G  Z,  a0  f  0,  such  that  the  representation  is  complete 
and  stable.  Of  particular  interest  is  the  sequence  of  scales  where  the  elementary 
dilation  step  a0  =  2  because  it  leads  to  an  octave  by  octave  partitioning  of  the 
frequency  domain. 
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In  order  to  reconstruct  s(t)  from  its  wavelet  transform  restricted  to  b  G  K  and 
a  =  a30,  j  G  Z.  the  Fourier  transform  of  must  satisfy  a  stricter  condition  than 
Equation  2.20  given  by  [Chu92] 

+00  2 

^<  Yl  K;(°o7)    <  B,  (2.22) 
j=— 00 

for  almost  all  7  G  R  for  some  constants  A  and  £  with  0  <  ,4  <  B  <  00.  Reconstruc- 
tion is  possible  by  using  the  following  formula  [Chu92] 

»(*)  =  E  /+°° {«o  W)K&)}  {«o>  (^r) }  db, 


where  the  Fourier  transform  of  r{t)  is  given  by 


•0(7)  = 


+00 


j=-oo 


If  satisfies  the  so-called  "stability  condition"  in  Equation  2.22  with  a0  =  2,  then 
ij}{t)  is  called  a  dyadic  wavelet  [Chu92]. 

In  addition  to  sampling  the  scale  parameter  a  =  aJ0,  j  G  Z,  it  is  possibl(>  to 
choose  a  sampling  of  the  translation  parameter  b  and  still  obtain  a  complete  and 
stable  representation.  Intuitively,  the  discretization  of  /;  has  to  be  chosen  such  that 
ipa;b{t)  "covers"  the  whole  time  axis  for  each  scale  a  =  aJ0,  j  G  Z.  From  Equation  2.15 
we  have  that  the  standard  deviation  of  ipa.b(t)  is  j«|  times  the  standard  deviation  of 
tp(t).  Hence,  in  order  for  ipa-,b(t)  to  "cover"  the  whole  time  axis  at  scale  a  =  a3Q,  the 
discretization  of  b  has  to  be  proportional  to  \a\.  It  follows  that  at  scale  a  =  aJ0  the 
translation  has  to  be  b  =  kb0aJ0,  where  bQ  >  0  is  the  elementary  translation  step.  This 
leads  to  a  discretized  family  of  wavelets  given  by 


r„(/)  =  a-fy  (  ]  =  flo  H  (aojt  -  kbQ) . 


"0 


In  order  to  reconstruct  s(t)  from  its  wavelet  transform  restricted  to  b  = 
kb0aJQ,  j,  k  G  Z  and  a  =  aJ0,  j  G  Z,  tp(t)  must  satisfy  a  stricter  condition  than  Equa- 
tion 2.22.    The  admissibility  condition  for  this  reconstruction  is  the  existence  of 
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.4  >  0,  B  <  oo  so  that 


A\\s\\22<       \WWiAt))\2  <  B\\a\\ 


2 
2 


for  all  s(t)  e  L'2{R)  [Dau92].  In  other  words,  the  family  of  functions  7/^(0  defined 
in  Equation  2.23  constitute  a  frame.  A  brief  review  of  frames  is  presented  in  Section 
2.2.1. 

Given  a  frame  of  wavelets,  reconstruction  of  s(t)  from  the  (s(t),ip^k{t))  is 
possible  by  using  the  following  formula 


where  the  family  of  functions  ^(t),  j,  A;  G  Z  is  called  the  dual  frame  ofipjtlc(t),  j,  k  e 
Z  (see  Section  2.2.1). 

If  the  dual  frame  of  7/>jifc(rj),  3,k  G  Z  is  of  the  form 


then  Equation  2.24  is  referred  to  as  a  wavelet  series  and  tp(t)  is  known  as  the  dual 
wavelet  of  ip(t).  Given  a  wavelet  frame,  the  existence  of  a  dual  wavelet  is  important  for 
computational  purposes.  Below,  we  overview  the  classification  of  wavelets  according 
to  orthogonality  and  highlight  the  existence  of  a  dual  wavelet  in  each  case.  A  wavelet 
can  be  classified  according  to  orthogonality  as  follows: 

1.  A  wavelet  ip(t)  is  called  semiorthogonal  if  it  satisfies  (tpj,k{t),^i,m{t))  =  0,  j  ^ 

I.  j,  k,  /,  777  G  Z. 

Remarks:  It  can  be  shown  [Chu92]  that  for  every  semiorthogonal  wavelet,  ip(t), 
there  exist  a  dual  wavelet,  xp(t),  such  that  the  pair  (ij)(t),^(t))  satisfies  the 


•s(/)  -  £   £  <*(/)•  r,,(0>^'(0. 


(2.24) 


where 
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biorthogonality  property 

(V>j,fc(<),  km{t))  =  Sjtl6k,m,j,  k,  I,  m  E  Z.  (2.25) 

2.  A  wavelet  ip(t)  is  called  nonorthogonal  if  it  is  not  a  semiorthogonal  wavelet. 

3.  A    semiorthogonal    wavelet    ip(t)    is    called    orthogonal    if    it  satisfies 

</W0)  =  SjjtSkjm,  h  k,  I,  m  €  Z. 

Remarks:  It  can  be  shown  that  an  orthogonal  wavelet  is  self  dual. 

4.  A  nonorthogonal  wavelet  ^{t)  is  called  biorthogonal  if  there  exist  a  dual  wavelet. 
i/>(i),  such  that  the  pair  (ip{t),ip(t))  satisfies  the  biorthogonality  property  in 
Equation  2.25. 

The  computation  of  Equation  2.24  can  be  a  burden  for  nonorthogonal  wavelets 

for  which  a  dual  wavelet  does  not  exist.  In  such  a  case,  it  is  advantageous  to  work  with 

frames  which  are  almost  tight  ("snug  frames"),  i.e.,  frames  which  have  B/A  - 1  «  1, 

because  the  ^k(t)  can  be  approximated  by  ipjtk(t)  (see  Section  2.2.1  and  [Dau92]  for 

more  details). 
2.2.1  Frames 

In  this  section  we  present  a  brief  review  of  frames.  For  a  more  detailed  and 
rigorous  account  of  frames  see  [Dau92] . 

A  family  of  functions  ipj,  j  G  J  in  a  Hilbert  space  U  is  called  a  frame  if  there 
exist  ,4  >  0.  B  <  oc  so  that  for  all  /  G  %, 

^||/||2<El^^)|2<5|l/H2,  (2-26) 

where  A  and  B  are  called  the  frame  bounds.  If  the  two  frame  bounds  are  equal,  the 
frame  is  called  a  tight  frame.  In  a  tight  frame  we  have,  for  all  /  e  H, 

£|</,¥>;>|2  =  4I/II2, 


1 1 


which  implies 

A{f,g)  =J^</,  Vi><Vi,P>, 
or  (at  least  in  the  weak  sense  [Dau92]), 

f  =  A-^(f^j)^r  (2-27) 

J'GJ 

Although  Equation  2.27  is  reminiscent  of  the  expansion  of  /  into  an  orthonormal 
basis,  it  is  important  to  note  that  frames,  even  tight  frames,  are  not  (orthonormal) 
bases.  The  family  of  functions  tpj:  j  e  J  are  typically  not  linearly  independent.  If 
the  frame  is  tight,  and  if  ||^||  =  1  for  all  j  €  Jf,  then  A  =  B  gives  the  "redundancy 
ratio."  If  this  ratio  equals  to  1,  then  the  tight  frame  is  an  orthonormal  basis  [Dau92]. 

Equation  2.27  gives  a  trivial  way  to  recover  /  from  the  (/,  ty),  if  the  frame  is 
tight.  Consider  now  recovering  /  from  frames  that  are  not  tight.  Let  us  define  the 
frame  operator  F  from  H  to  £2(JT)  as 

Since  ipj  constitute  a  frame,  it  follows  from  Equation  2.26  that  ||F/||2  <  5||/||2,  that 
is,  F  is  bounded,  which  means  it  is  possible  to  find  its  adjoint  operator  F*  [Dau92]. 
The  adjoint  operator  F*  of  F  can  be  computed  from  the  following  relation 

(F*cJ)   =   (c,Ff)=  $><7^> 

so  that 

at  least  in  the  weak  sense  [Dau92].  From  this  it  follows  that 

E  </'  V>i)  Vi  =  F*Ff. 
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Thus,  the  definition  of  F  implies 

£l(/^>l2Him2H^^/>- 

Hence,  in  terms  of  F,  the  frame  condition  in  Equation  2.26  can  be  written  as 

AI<F*F<BI.  (2.28) 

where  /  is  the  identity  operator  [Dau92].  This  implies,  in  particular,  that  F*F  is 
invertible  (see  Lemma  3.2.2  in  [Dau92]),  and  that  the  operator  (F*F)_1  satisfies, 

B~l  I  <  (F'F)"1  <  A~lI. 

Applying  the  operator  (F*F)_1  to  the  family  of  functions  ipj,  j  G  J,  leads  to 
another  family  of  functions  denoted  by  tf>j,  j  G  J  where 

<Pj  =  (F*F)"Vi,  3  G  J. 

The  family  of  functions  ipj,  j  G  Jf  also  constitutes  a  frame  with  frame  bounds  B~l 
and  A~x  [Dau92],  that  is 

B-1||/||a<Sl</,^>|2<A-1||/||a. 

It  can  be  verified  (see  Proposition  3.2.3  in  [Dau92])  that  the  associated  frame  operator 
F  from  H  to  p(J)t  (Ff)j  =  satisfies  F  =  F{F*F)~l,  F*F  =  (F*F)-\ 

F*F  =  I  =  F*F  and  FF*  =  FF*  is  the  orthogonal  projection  operator,  in  ^2(Jf), 
onto  Ran(F)  =  Ran(F).  The  family  of  functions  ipj,  j  G  J  is  called  the  dual  frame 
of  ipj,  j  G  J.  It  is  easy  to  verify  that  the  dual  frame  of  <pj,  j  G  J  is  (fj,  j  G  Jf.  From 
F*F  =  I  -  F*F  it  follows  that 

£  </,^>  ft  =  /  =  £</,  (2-29) 

Hence,  we  have  a  way  to  recover  /  from  the  (/,  tpj),  where  the  only  thing  we  need 
to  do  is  to  compute  the  (pj  =  (F*F)~1(pj.   Note  that,  in  general,  if  the  frame  is 
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redundant,  there  exist  other  functions  in  W.  that  could  equally  well  play  the  role  of 
the  ipj  and  lead  to  a  reconstruction  formula.  This  follows  from  the  fact  that  the 
ipj  are  not  linearly  independent  in  the  general  case.  Equation  2.29,  however,  yields 
the  most  "economical"  representation  of  /  in  the  following  sense  (see  Proposition 
3.2.4  in  [Dau92]):  Consider  the  first  half  of  Equation  2.29,  /  =  52ja(f,<pj)<pj.  If 
/  =  YLjes  (/'  Vj)  u3  f°r  some  other  family  of  functions  Uj,  j  G  J.  then  it  can  be  shown 
[Dau92]  that  £V€j  |(/,Uj)|2  >  £i€,  \(f><Pj)\2-  Similarly,  consider  the  second  half  of 
Equation  2.29,  /  =  £ -eJI  (f,<Pj)  <Pj-  K  ./'  =  Y^jeJcjVj  f°r  some  c  e  ^2U)'  an(I  ^  not 
all  cj  equal  (f,<Pj),  then  it  can  be  shown  [Dau92]  that  J2jel  lci|2  —  ^ZjeJ  K/»^i)l  ■ 

Computing  the  (pj  involves  the  inversion  of  F*F.  If  B  is  close  to  .4  (r  = 
B/A  -  1  «  1),  then  from  Equation  2.28  we  have  that  F*F  is  "close"  to  ±&  I.  This 
implies  that  (F*F)_1  is  "close"  to  j^gl,  and  that  is  "close"  to  -j^gVi  [Dau92]. 
This  motivates  the  following  reconstruction  formula  for  /  [Dau92]. 

/  =  T^B  S (/'  ^>  ^  +  R^  (2:30) 

j€l 

where 

R  =  /  -  — (2-31) 

It  follows  from  Equations  2.28  and  2.31  that  -f^£/  <  i2  <  f^f  /.  This  implies 
that  \\R\\  <  ^|  =  2+7  <  1-  If  ''  is  small,  the  rest  term  Rf  in  Equation  2.30  can 
be  dropped,  leading  to  a  reconstruction  formula  for  /  which  is  accurate  up  to  an 
L2-error  of  j^H/H  [Dau92].  Even  if  r  is  not  small,  it  is  possible  to  write  an  algorithm 
for  the  reconstruction  of  /  with  exponential  convergence  [Dau92].  From  Equation 
2.31  we  have  that 

.4  +  B 

This  implies  that 

{rF)~1  =  A-h{I-RrK 
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Since  \\R\\  <  1,  the  series  YlT=o^k  converges  in  norm,  and  its  limit  is  (/  -  R)  1 
[Dau92].  It  follows  that 


<Pj  =  (F*F)-Vi 


A  +  B 


(2.32) 


Note  that  the  zeroth  order  term  in  the  above  equation  leads  exactly  to  Equation  2.30 
with  the  rest  term  dropped.  Better  approximations  can  be  obtained  by  truncating 
after  N  terms  [Dau92], 

2 


v 


A  +  B 


k=0 


with 


/  -  E  </'  vi>  $ 


< 


2  +  r 


N  +  l 


which  becomes  exponentially  small  as  N  increases,  since  ^  <  1.  In  particular,  the 
ifj  can  be  computed  by  an  iterative  algorithm  [Dau92], 


Similarly,  /  can  be  computed  by  an  iterative  algorithm  [Dau92], 


with 


f  =  (F*F)-1(F*F)f=  lim 

iv— >oo 


9 

/.v  =  /jv-i  +  -J^Tb  5^  ^'  ^  ~  ^N~l,lfj^ 


2.2.2    Wavelet  Frames 


j€3 


Not  all  choices  for  tp(t),  aQ,  b0  lead  to  frames  of  wavelets,  even  if  ip(t)  is  admis- 
sible. Daubechies  [Dau92]  gives  some  general  conditions  on  ijj(t),  a0,  b0  under  which 
a  frame  is  obtained  and  derives  estimates  for  the  frame  bounds.  These  results  can  be 
summarized  as  follows  (see  Proposition  3.3.2  in  [Dau92]):  If  ip{t),  a0,  are  such  that 


inf     E  h/;(ao7) 

l<|7|<a0  I 
j  =  -oo 

sup  E  mao"f) 

1<|7|<«0  1 


>o, 


<  oo, 


(2.33) 


is 


and  if 


/j(7s)  =    sup     V  ^(oo7) 

l<|7|<a0  " 


t/>(aJo7  +  Is) 


]=-<x> 


decays  as  fast  as  (1  +  |7s|)"(1+f).  with  e  >  0.  then  there  exist  blln  >  0  such  that 
the  (.-/,(/)  constitute  a  frame  for  all  choices  b0  <  bihr.  For  b0  <  bthr,  the  following 
expressions  are  frame  bounds  for  the  ipj,k(t) 


+00  r 


A  =  r\j*l  E  I ■  E 


fro   1  K|7l<ao 

j  =  -oc 


A"  =  —  oo 
A'  Ml 


1  1 


n  = 


y\  sup  e        +E  nrk M 


2tt 


i  i 


The  condition  on  /?(7s)  and  Equation  2.33  are  satisfied  if  |V>(7)|  <  Cb'|Q(l  +  \"f\)~0 
with  a  >  0,  0  >  a  +  1.  It  can  also  be  verified  that  if  the  <.>(/)  in  Equation  2.23 
constitute  a  frame  with  frame  bounds  A  and  B  for  some  choice  of  a0,  b0,  then  the 
Fourier  transform  of  i>{t)  satisfies 


m  <  E  r(a°7) 


<  bnB 


(2.34) 


almost  everywhere  [Chu92]. 

As  mentioned  before  it  is  advantageous  to  set  the  elementary  dilation  step 
o0  =  2  because  it  leads  to  an  octave  by  octave  partitioning  of  the  frequency  domain. 
Hereafter,  we  will  assume  this  partitioning  of  the  frequency  domain  unless  stated 
otherwise. 

2.2.3    Multivoice  Wavelet  Frames 

In  Section  2.2  we  mentioned  that  it  is  advantageous  to  use  wavelet  frames 
with  B/i-Kl,  because  the  ^k(t)  can  be  approximated  by  ipjjt(t).  This  implies 
that  the  sum  in  Equation  2.34  should  be  almost  constant  for  7  ^  0.  This  imposes 
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a  strong  restriction  on  ip(t)  that  is  not  generally  satisfied.  In  order  to  overcome  this 
problem  without  having  to  give  up  too  much  freedom  in  choosing  ij)(t)  or  its  band- 
width, one  can  adopt  the  use  of  different  voices  per  octave  (i.e.,  suboctave  sampling) 
[Dau92].  This  can  be  achieved  by  using  different  wavelets  tpi{t),  ip2(t),  ■■■  ,*l>ft(t)  and 
looking  at  the  frame  tpv.j,k(t)  =  -ijtv  -  fc&o)  ,  j >*  €  Z.  v  =  1, . . .  ,  .V.  In  [Dau92], 
Daubechies  gives  the  following  expressions  for  the  frame  bounds  of  this  multivoice 
frame 


C  N        +  00  V  r 

u    I,  v=l  j=-oo  fc#0  i»=l  L 


1  1/2' 


(  JV      +00  AT     r  , 

r  ™p  E  E  i^7)i2  +  £E  a( 


2tt 


-i  1/2 


with 


+  00 


&(7s)  =   sup    ^  |^(2J7)||^(2J7  +  7.)I- 


i<H<2J=_00 


It  is  suggested  in  [Dau92]  that  by  choosing  the  $1(7),...  ,^n{i)  to  have  slightly 
staggered  frequency  localization  centers,  coupled  with  good  decay  at  00,  one  can 
achieve  B/A  -  1  <C  1.  One  choice  favored  by  several  authors  [Gro90,  Dau92.  Lee96] 
is  to  take  "fractionallv"  dilated  versions  of  a  single  wavelet 


Mt)  =  2 


(2.35) 


In  this  case 


N  +00 


+00 


|2  (2.36) 

r=i  j=— 00  j=— 00 

which  can  be  made  almost  constant  by  choosing  JV  large  enough  [Dau92]. 


CHAPTER  3 
DISCRETE  FILTERS  AND  WAVELETS 

In  this  chapter  we  review  the  computation  of  the  wavelet  transform  and  its 
inverse  for  nonorthogonal  wavelets.  Specifically,  we  review  the  computation  of  the 
wavelet  coefficients  in  equation  2.24  by  using  discrete  filters  and  the  inverse  transform 
in  a  nonorthogonal  multiresolution  framework.  We  will  refer  to  equation  2.24  loosely 
as  a  wavelet  "series"  even  in  the  case  that  the  dual  frame  is  not  generated  by  the 
dilations  and  translations  of  a  single  function.  We  will  also  use  the  term  nonorthogo- 
nal wavelet  loosely  to  refer  to  a  nonorthogonal  wavelet  for  which  a  wavelet  dual  does 
not  exist.  We  also  restrict  the  scope  of  our  review  to  the  case  where  the  elementary 
dilation  step  a0  =  2  and  the  elementary  translation  step  b0  =  1. 

3.1    The  A  Trous  and  Mallat's  Algorithms 

Recall  from  Section  2.2  that  the  wavelet  transform  of  a  signal  s(t)  is  given  by 

/+00  /  ±  _  f  \ 

The  wavelet  coefficients  {s(t),ipj>k(t))  in  Equation  2.24  for  a0  =  2,  b0  =  1  can  be 
obtained  by  sampling  the  wavelet  transform  of  s(t)  with  a  =  2J .  j  G  Z  and  b  = 
2Jk,  j,k  G  Z.  We  will  refer  to  these  coefficients  as  the  decimated  wavelet  "series" 
coefficients  and  will  denote  them  by 

(Cts)hk  =  (VM(2*,2*fc) 

=   (s(t)^hk(t)) ,  (3.1) 

where  ipj,k{t)  =  (57  -  k)  is  the  generator  of  a  family  of  decimated  wavelets. 
We  will  refer  to  the  operation  of  computing  the  above  wavelet  coefficients  as  the 
decimated  wavelet  transform. 
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For  many  applications  the  above  sampling  is  inadequate.  In  particular,  a 
finer  grid  is  needed  when  computing  the  wavelet  transform  for  applications  that 
require  "good"  time-frequency  localization.  This  can  be  achieved  by  computing  an 
undecimated  wavelet  transform  and  the  addition  of  voices  [She92].  The  undecimated 
wavelet  transform  of  a  signal  s(t)  is  obtained  by  sampling  the  wavelet  transform  of 
s(t)  with  a  =  2j,  j  6  Z  and  =  k,  k  E  Z.  We  will  refer  to  these  coefficients  as  the 
undecimated  wavelet  "series"  coefficients  and  will  denote  them  by 

(Mi,*   =  (VM(2M;) 

=   (s(t).ii(t)),  (3.2) 

where  ip{(t)  =  -j^ip  (^)  is  the  generator  of  a  family  of  undecimated  wavelets. 

In  what  follows  we  will  review  two  separately  motivated  implementations  of 
the  wavelet  transform,  the  algorithme  a  trous  [Hol90]  and  Mallat's  [Mal89a]  multires- 
olution  decomposition.  These  algorithms  are  both  special  cases  of  a  single  filter  bank 
structure,  the  discrete  wavelet  transform  [She92],  the  behavior  of  which  is  governed 
by  one's  choice  of  filters.  The  a  trous  algorithm  was  originally  devised  as  a  com- 
putationally efficient  implementation  of  the  wavelet  transform;  however,  it  is  more 
properly  viewed  as  a  nonorthogonal  multiresolution  decomposition  for  which  the  dis- 
crete wavelet  transform  computes  a  "sampled"  continuous  wavelet  transform  (wavelet 
"series"  coefficients)  exactly  [She92].  On  the  other  hand,  Mallat's  algorithm  was  de- 
vised as  an  orthogonal  multiresolution  decomposition  [Mal89a]  that  when  initialized 
properly  computes  a  "sampled"  continuous  wavelet  transform  [She92].  For  a  more 
detailed  and  rigorous  account  of  these  algorithms  see  [Mal89a.  Hol90,  She92]. 

The  following  additional  terminology  and  notation  will  be  used  throughout 
this  chapter.  The  discretized  wavelet  transform  (with  sampling  period  T  =  1)  is 
given  by 

{w^)(aM  =  ^XS[n)Vi^)-  (3'3) 
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The  discretized  decimated  wavelet  "series"  coefficients  are  given  by 

M)M  =  KS)(2'>'fc)     ^.sWrg  -  k).  (3.4) 
Similarly,  the  discretized  undecimated  wavelet  "series"  coefficients  are  given  by 

Mn  =  K*)(2J>  *)  =  ^=  E  s(n)  *  (^)  ■  (3-5) 

Signals  and  filters  in  bold  face  will  be  treated  as  vectors.  The  A;th  element 
of  vector  /  is  denoted  by  [f]k  =  fk.  The  symbol  t  will  be  used  for  the  mirror  filter 
[/tjfc  _  yt  _  J  ^  Tne  decimator  operator  is  denoted  by  a  matrix 

Dk,m   =  S{2k  -  m) 

where  8kjn  is  the  Kronecker  delta  and  [6}k  =  6(k)  =  Skt0.  The  dilation  operator  is 
denoted  by  a  matrix 

Uk,m  =  8(k  -  2m) 

Convolution  of  /  and  s  is  denoted  by  [/  *  s]k  =  Y,m  fk-mSm-  Convolution 
followed  by  decimation  becomes  [£>(/  *  s)]k  =  Y.mDk,m[f  *  s]m  =  [/  *  shk  = 
]£m  /2fc-m«m-  The  following  shorthand  will  be  used  for  the  convolution  of  /  and  s, 
f*s  =  Fs,  where  Fm,„  =  fm-„.  It  is  easy  to  verify  that  D{f*s)  =  DFs.  The  symbol 
f  will  also  be  used  for  the  Hermitian  transpose  of  matrices  [J1*]™^  =  / 

n—m  fm—w 

For  discrete  signals  the  translation  operator  is  denoted  as 

[7n,s]n  —  .S'm-ni 

and  for  continuous  signals  as 

(TTs)(t)  =  s(t-r). 
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The  z  transform  of  a  discrete  signal  sn  is  given  by 

n 

Finally,  we  state  the  following  two  identities  [She92],  which  will  be  used  later.  For 
any  F  of  the  form  Fm,„  =  fm-n  we  have  that 

[(DFYU  =  [(DF)%^2Jn  (3.6) 

and 

^[(DFy]n,kzk=z»nl[f(z2r).  (3.7) 
3.1.1    The  Decimated  A  Trous  Algorithm 

The  a  trous  algorithm  was  originally  devised  as  an  efficient  implementation  for 
computing  wavelet  series  coefficients.  To  achieve  this,  the  wavelet  series  coefficients  in 
Equations  3.1  and  3.2  are  approximated  by  their  discretized  counterparts  in  Equations 
3.4  and  3.5,  respectively.  As  a  starting  point  consider  implementing  Equation  3.4. 
Clearly,  as  j  increases  ip(t)  must  be  sampled  at  progressively  more  points,  creating  a 
large  computational  burden  [Hol90,  She92].  The  solution  proposed  by  [Hol90]  is  to 
approximate  nonintegral  points  via  a  finite  filter  f.  As  an  example,  let  f]  be  the 
filter  ^=(0.5,1.0,0.5).  Then, 

/^V-  ,t  _    /  0(i)>  for  n  even, 

V2      /»-»Vl*J   "    {  I  jy(n=i)  +  ^(»±1))  ,   for  n  odd, 


approximates  a  sampling  of  tp{t/2).  Let  g  be  a  filter  defined  by  gn  =  i'(-n).  Then 
the  above  interpolation  can  be  computed  by  first  dilating  g]  and  then  convolving  the 
result  with  a  filter  /f  which  leaves  the  even  points  fixed  and  interpolates  to  get  the 
odd  points  [She92].  The  condition  that  /  be  the  identity  on  even  points  is  sufficiently 
important  that  merits  the  following  definition:  A  lowpass  filter  /  is  said  to  be  an  a 
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trous  filter  if  it  satisfies 

f2k  =  8(k)/V2.  (3.8) 

The  y/2  is  simply  a  convenient  means  of  including  the  normalization  factor  of  Equa- 
tion 3.4  in  the  filter. 

It  can  be  verified  [She92]  that  the  interpolation  operation  can  be  written  as 

follows 

\f*(Ug%   =  [F^UgX 

=  fn-2k9k 
k 

k 

This  result  and  Equation  3.4  leads  to  the  following  approximation  for  the  discretized 
decimated  wavelet  series  coefficients  at  j  =  1, 


Mhk  =  -7|5Z8(n)^(f -*) 

*  n 

~  sn         f  n-2k-2m9m 

n  m 

—  Qk—m        fim— n^n 

m  n 

=   [g*  D{f  *  s)]fc 
=  [GDFs]k, 

where  [s]n  =  sn  =  s(n).  Proceeding  inductively,  one  can  find  the  following  approxi- 
mation for  the  discretized  decimated  wavelet  series  coefficients  for  all  j  [She92] 

(c*a)it*  «  (3.9) 
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I  2 


Sj+i 


Figure  3.1:  Block  diagram  of  one  stage  of  the  discrete  decimated  wavelet  transform. 

Remarkably,  the  right-hand  side  in  the  above  equation  is  the  discrete  decimated 
wavelet  transform  of  s  at  level  j  [She92],  which  we  denote  as  follows1 


[Cj(8)]k  =  [G(DFYs},. 


(3.10) 


Equations  3.9  and  3.10  yield  the  following  formula  for  the  discretized  decimated 
wavelet  series  coefficients 


(3.11) 


The  discrete  decimated  wavelet  transform  can  be  computed  iteratively  over  j  as 
follows 

[8j+i]k     =     /]  /2fc-n[gj]n 
n 

=  [D(f*8j)]k 


n 


(3.12) 


where  [s0}k  =  sk.  A  block  diagram  of  one  stage  of  the  transform  is  shown  in  Figure 
3.1.  Except  for  decimation  of  the  output  (see  Section  3.1.3)  the  a  trous  algorithm 
described  in  [HolOO]  is  given  by  Equations  3.11  and  3.12. 

The  original  a  trous  algorithm  made  no  statements  regarding  the  accuracy  of 
the  approximation  in  Equation  3.9  or  even  the  discretization  from  Equation  2.9  to 


'The  notation  [c,-]*  =  c,.*  will  be  used  as  a  shorthand  for  [Cj(a)]fe. 
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Equation  3.3.  A  major  step  forward  towards  treating  this  question  lies  in  the  results 

of  [She92]  outlined  in  Section  3.2. 
3.1.2    Mallat's  Algorithm 

In  this  section  we  review  Mallat's  multiresolution  algorithm.  For  a  more  de- 
tailed and  rigorous  account  of  this  algorithm  see  [Mal89a,  Vet92,  She92].  Mallat's 
algorithm  has  basically  the  same  structure  as  Equations  3.12,  that  is, 


where  [s0]k  =  sk.  In  keeping  with  the  literature,  the  lowpass  filter  is  denoted  by  h 
instead  of  /.  The  constraints  on  h  and  g  which  ensure  an  orthonormal  multiresolution 
analysis  are  [Mal89a,  She92]: 

1.  Perfect  reconstruction 


[D(h  *  8j)]k 


[D(g*8j)]k 


(3.13) 


(3.14) 


2.  Orthogonality  of  h  and  g  with  respect  to  even  shifts 


(3.15) 


3.  Bandpass  condition  on  g 


(3.16) 


n 


4.  Lowpass  condition  on  h 


(3.17) 
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In  matrix  notation  Equations  3.14  and  3.15  may  be  written  as  follows 

(WU){DH)  +  {G]U){DG)  =  I.  (3.18) 

(DH)(G*U)  =  0,  (3.19) 

where  Hnhn  =  hm-n.  Multiplying  Equation  3.18  on  the  left  by  DH  and  using 
Equation  3.19  we  have  that 

(DH){H^U)  =1. 

It  follows  that  h  and  its  shifted  versions  by  even  shifts  form  an  orthonormal  set 

^  h2m-jh2n-j  =  6m,n,  m,  fl  <E  Z. 
j 

Similarly,  multiplying  Equation  3.18  on  the  left  by  DG  and  using  Equation  3.19  we 
have  that 

(DG){GW)  =  I 

It  follows  that  g  and  its  shifted  versions  by  even  shifts  also  form  an  orthonormal  set 

?J  92m-j92n-i  =  <Wi,  TH,  U  G  Z. 
J 

These  two  results  together  with  Equation  3.15  imply  that  Equations  3.13  is  an  or- 
thogonal decomposition  of  the  discrete  signal  sr  Furthermore,  Equation  3.14  implies 
that 

[8j)k  =  i(&U)sJ+i  +  (&U)dj+l]k. 

From  Equations  3.14  3.17  it  follows  that  Equations  3.13  represents  a  wavelet 
decomposition  as  described  below  [Mal89a,  She92]. 

Define  a  scaling  function  <p(t)  with  Fourier  transform  given  by 

r=l  V  ~ 
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It  follows  that 

which  in  the  time  domain  takes  the  form 

cj)(t)  =  ^2~h_kV2(f)(2t  -  k). 
k 

Therefore,  the  dilates  and  translates  of  </>(£), 
have  the  property  [She92] 

n 

Define  a  wavelet  function  by 

il>(t)  =  Y,g_k>/2<i){2t-k).  (3.20) 

Then,  using  Equation  3.15  and  the  above  properties,  one  can  show  that  the  family 
of  wavelets 

are  orthonormal,  and  that  the  dj  are  the  coefficients  of  the  expansion  of  s(t)  in  terms 
of  i[>j,k{t),  that  is 

provided  that 

/ +  00   
s(i)<£(t  -  fc)<ft. 
-00 

Therefore,  Mallat's  algorithm  when  initialized  properly  (see  the  equation  above)  com- 
putes exactly  a  "sampled"  continuous  wavelet  transform  [She92]. 
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3.1.3    The  Undecimated  A  Trous  Algorithm 

As  mentioned  before  for  many  applications  the  time-frequency  sampling  grid 
provided  by  the  decimated  wavelet  transform  is  inadequate  when  "good"  time- 
frequency  localization  is  required.  One  way  to  overcome  this  limitation  is  by  com- 
puting an  undecimated  wavelet  transform.  Following  the  same  approach  described 
above  the  wavelet  series  coefficients  in  Equation  3.2  can  be  approximated  by  their 
discretized  counterpart  in  Equation  3.5.  From  Equations  3.4  and  3.5  it  follows  that 
(cv,s)M  =  (utys)(2J',2J'A;)  and  (i>s)j,fe  =  (uty«)(2>,fc)  should  coincide  at  k  =  0.  Uti- 
lizing this  fact,  one  can  obtain  the  kth  undecimated  wavelet  series  coefficient  (i>s)jifc 
by  translating  the  signal  back  by  k  and  computing  the  discretized  decimated  wavelet 
series  coefficient  (c*(T_fcs))i|0.  This  result  together  with  Equation  3.11  yields  the 
following  approximation  for  computing  the  discretized  undecimated  wavelet  series 
coefficients 

{rts)j,k   =  MT-ks))jfi 

«  [c,(T_*a)]0l  (3-21) 

where  [s}n  =  sn  =  s(n).  Remarkably,  the  right-hand  side  in  the  above  equation  is 
the  discrete  undecimated  wavelet  transform  of  s  at  level  j  [She92],  which  we  denote 
as  follows2 

[rj(8)]k  =  [cJ(r_Jt*)]0.  (3-22) 

Equations  3.21  and  3.22  yield  the  following  formula  for  the  discretized  undecimated 
wavelet  series  coefficients 

Mw  «  Ma)]*.  (3-23) 
2The  notation  [rj]k  =  rj,k  will  be  used  as  a  shorthand  for  [rj(s)]k 
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A  few  observations  about  the  decimated  and  undecimated  discrete  wavelet 
transforms  are  in  order.  First  notice  that  in  general  c,  is  not  translation  invariant 

[c,(Tma)]fc  =  ^[G(DFy']M[i-m 

=  Y,iG(DFykn+m[s}n 

n 

/  Y^[G(DF)j}k-m,n[a]n. 

However,  if  one  replaces  m  by  2jm  in  the  above  equation  and  uses  Equation  3.6.  the 
last  step  becomes  an  equality  [She92] 

[Cf(T*m«)]*   =   M»)]*-m-  (3-24) 

Thus,  translating  s  by  2jm,  translates  octave    by  m. 

On  the  other  hand,  it  can  be  verified  that  r3  is  translation  invariant 

[rj(Tms)]k  =  [c,-(Tm_fc*)]o 
=  [rj(s)]fc_m. 

Also,  from  Equations  3.22  and  3.24  it,  follows  that  sampling  rjtk  every  23  points 
produces  exactly  Cjik,  that  is, 

Cj,fc  =  fj,2ik- 

An  iterative  algorithm  for  computing  the  discrete  undecimated  wavelet  transform  can 
be  obtained  by  taking  z  transforms  [She92].  From  Equations  3.10  and  3.22 

X>,.,c  k   =  ££[G(DF)']0,,nW2-* 

k  m  k 

in 

Applying  Equation  3.7  to  the  above  equation  it  follows  that 
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Ujg 


Figure  3.2:  Block  diagram  of  one  stage  of  the  discrete  undecimated  wavelet  transform. 

where  j  =  0  is  understood  to  mean  there  are  no  factors  of  f(z). 

It  is  easy  to  see  that  Uj  f  is  /  with  2J  -  1  zeros  inserted  between  every 
pair  of  filter  coefficients  and  that  its  z  transform  is  f{zp).  Therefore  the  discrete 
undecimated  wavelet  transform  can  be  computed  iteratively  over  j  as  follows 

[sj+i]k   =  [(Uif)*Sj)]k 

N*   =  [(Uig)*aj]k 

where  [sQ]k  =  sk.  The  above  result  together  with  Equation  3.23  is  essentially  the  a 
trous  algorithm  described  in  [HolQO].  A  block  diagram  of  one  stage  of  the  transform 
is  shown  in  Figure  3.2. 

3.2    The  A  Trous  Algorithm  as  an  Exact  Wavelet  Transform 

In  the  previous  section  we  analyzed  two  separately  motivated  implementations 
of  the  wavelet  transform.  We  learned  that  the  a  trous  algorithm  was  originally  devised 
as  a  computationally  efficient  approximation  of  the  "sampled"  wavelet  transform 
(Cil>s)j,k  whose  implementation  can  be  viewed  as  a  discrete  wavelet  transform  for 
which  the  filter  /  is  an  interpolator  satisfying  Equation  3.8  and  the  filter  g  is  obtained 
by  sampling  the  wavelet  function  ij;(t).  We  also  learned  that  Mallat's  algorithm 
is  a  discrete  wavelet  transform  with  filters  h  and  g  subject  to  the  constraints  in 
Equations  3.14  3.17  that  when  initialized  properly  computes  exactly  a  "sampled" 
continuous  wavelet  transform  for  the  wavelet  function  i)(t)  in  Equation  3.20.  Next 
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we  will  review  the  results  of  [She92]  which  show  that  the  a  trous  algorithm  can 
be  viewed  as  a  nonorthogonal  multiresolution  decomposition  for  which  the  discrete 
wavelet  transform  computes  a  "sampled"  continuous  wavelet  transform  exactly. 
Define  a  scaling  function  <j>(t)  with  Fourier  transform  given  by 

OO 

r=]  V- 

or  equivalent!}' 

</>(t)  =  lira  y\{DF)%,kV2iX{2H  -  k),  (3.26) 
k 

where 

y(t)  =  S  1   forte  [-1/2,1/2), 
\  0  otherwise. 

For  normalization  purposes  f{eiu3)\w=Q  =  \/2,  which  implies 

£  /,  =  s/2  (3.27) 

k 

and 

-+oo 

<j)(t)dt  =  0(0)  =  1. 


/' 

j  — < 


Suitable  regularity  conditions  for  the  inverse  Fourier  transform  of  the  product  in 
Equation  3.25  to  converge  to  a  reasonably  behaved  function  may  be  found  in  [Dau88, 
Mal89b,  She92].  One  can  verify  [She92]  that  the  dilates  and  translates  of  (j>(t)  have 
the  property 

<f>j+hk(t)   =  53[SF)fc,„^,«(*). 

n 

=     £72*-n<M*)-  (3-28) 
n 

Furthermore,  if  /  is  a  finite  filter,  then  </>(£)  has  finite  support  [She92]. 
Define  the  wavelet  function 

=  Y,9kHt  +  *)  =  £  9k<t>(t  -  k).  (3.29) 
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We  will  refer  to  ij'{t)  as  the  approximation  of  ijj(t).  The  recursion  in  Equation  3.28 
implies  that 

<M*)  =  £[(^)j']m<M*). 

n 

This  expression  along  with  Equation  3.29  yields 

^(t)  =  ^[G(Df)j]M^(<). 

n 

Using  the  above  expression  in  Equation  3.4  we  have  that  the  discretized  decimated 
wavelet  "series"  coefficients  take  the  form 

=  £  «(0  £[g(W)']m<mo 

/  n 

n  I 

=  ^[G(DF)^.n^.S(l)^). 

n  / 

Define 

[a].  =  5>(05Fr^.  (3-30) 

then  the  above  equation  becomes 

{cir8)JJk  =  [G{DFYa]ki 

which  is  the  discrete  decimated  wavelet  transform  of  the  discrete  signal  sn  defined 
in  Equation  3.30.  Therefore,  using  the  initialization  in  Equation  3.30  the  discrete 
decimated  wavelet  transform  computes  exactly  the  discretized  decimated  wavelet 
"series"  coefficients  (c^s)^.    Similarly,  using  the  above  expression  for  i^k(t)  in 
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Equation  3.1  we  have  that  the  decimated  wavelet  "series"  coefficients  take  the  form 

/+00   
s(t)^k(t)dt 
■00 

=   j+0°  s(t)  \y^[G(DF)j}k,n<f>o,n(t)j  dt 


/+00   
s(t)<t>0,n(t)dt 
■00 
*  +  00 


/+oo   
s(t)<f>(t  -  n)dt. 

n 

Using 

/+O0   
s(t)<f>(t  -  n)dt  (3.31) 
-00 

the  above  equation  becomes 

(C^s)hk  =  [G(DFys}k, 

which  is  the  discrete  decimated  wavelet  transform  of  the  discrete  signal  sn  defined 
in  Equation  3.31.  Therefore,  using  the  initialization  in  Equation  3.31  the  discrete 
decimated  wavelet  transform  computes  exactly  the  decimated  wavelet  "series"  coef- 
ficients (Qs)^.  Notice  that  Equation  3.30  is  the  discretized  version  of  Equation 
3.31. 

Finally,  we  review  the  significance  of  the  a  trous  condition  in  Equation  3.8. 
One  may  show  that  a  filter  /  is  an  a  trous  filter  iff  0(n)  =  6(n)  [She92].  Using 
this  observation  it  can  be  shown  that  if  /  is  a  trous  then  (1)  Equation  3.30  becomes 
[«]„  =  s(n)  and  the  discrete  decimated  wavelet  transform  of  the  sampled  signal  s{n) 
yields  the  exact  discretized  decimated  wavelet  "series"  coefficients  (c^s)^  and  (2) 
the  wavelet  ip'(t)  defined  in  Equation  3.29  satisfies  ip' (ri)  =  g]n  [She92]. 

In  summary,  given  /  satisfying  Equations  3.8  and  3.27,  g  defined  by  g\  =  ij}(n) 
(where  tp(t)  is  an  arbitrary  wavelet),  and  given  the  initialization  in  Equation  3.31,  the 
discrete  wavelet  transform  computes  an  exact  "sampled"  continuous  wavelet  trans- 
form of  s(t)  using  the  wavelet  function  ip' (t)  defined  in  Equation  3.29.  Furthermore, 
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if  there  is  sufficient  regularity,  ip'(t)  and  if)(t)  will  be  ''close"  since  they  coincide  on 
the  integers  up  to  the  length  of  g. 

3.3    The  Undecimated  A  Trous  Algorithm  and  Voices 

As  mentioned  before  we  are  interested  in  computing  a  multivoice  undecimated 
wavelet  transform.  In  this  section  we  show  how  to  extend  the  algorithm  presented 
in  Section  3.1.3  into  the  multivoice  framework.  Recall  from  Section  2.2.3  that  in 
the  multivoice  case  instead  of  having  a  basic  wavelet  t/>(t)  one  deals  with  a  family  of 
wavelets  ipv(t),  v  =  1,...  ,N.  In  the  spirit  of  Equations  3.1  and  3.2  we  denote  the 
multivoice  decimated  and  undecimated  wavelet  "series"  coefficients  as 

(Qt,3)j)fc  =  (Wv,e5)(2^2^) 

and 

(*+s)jJt  =  mvs)(2j,k), 

respectively.  Similarly,  the  discretized  multivoice  decimated  and  undecimated  wavelet 
"series"  coefficients  are  given  by 


y  '  V2J 

and 


1  /  x  (n-k\ 

respectively.  Following  the  approach  presented  in  Section  3.1.1  let  gv  be  the  filter 
defined  by  gv;n  =  ipv{-n)  and  let  /  be  an  a  trous  filter.  Then 

~  [Cv;j(8)]k  (3-32) 

where 

[c*(a)J*  =  [Gv(DFys]k.  (3.33) 
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Notice  that  Equation  3.33  resembles  the  structure  of  a  discrete  wavelet  transform 
with  the  exception  that  we  now  have  a  family  of  filters  gv.  The  coefficients  [cv.j]k  = 
[Cy.j(s)]k  can  be  computed  iteratively  over  j  as  follows 

=  /j2k-n[Sj]n 
n 

=  [D(f*8j)]k 

[Cv;j]k     =  ^29v;k-n[Sj]n 

n 

=     [9v  *8j]k 

where  [s0]k  =  sk. 

Following  the  approach  presented  in  Section  3.2  let       be  defined  as  in  Equa- 
tion 3.25  and 

k  k 

Using  the  initialization  given  in  Equation  3.30  one  can  show  that  the  coefficients  on 
the  right-hand  side  of  Equation  3.32  are  the  exact  discretized  multivoice  decimated 
wavelet  "series"  coefficients  (c^s)j>  Furthermore  if  /  is  a  trous  the  result  [s]n  =  s(n) 
still  holds.  Similarly,  using  the  initialization  given  in  Equation  3.31  one  can  show 
that  the  coefficients  on  the  right-hand  side  of  Equation  3.32  are  the  exact  multivoice 
decimated  wavelet  "series"  coefficients  {C^  s)hk.  Again  if  /  is  a  trous  the  result 
ip'v{n)  =  ^v{n)  still  holds. 

Extending  the  above  results  to  the  undecimated  case  can  be  easily  accom- 
plished by  using  Equation  3.22,  that  is 

[rV:j{s)}k  =  [ct);J(r_fcs)]0, 

where  the  coefficients  [rVij]k  =  [rv.j(s)]k  can  be  computed  iteratively  over  j  as  follows 

[BJ+l]k     =  [(Ujf)*Sj)]k 
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Figure  3.3:  Block  diagram  of  one  stage  of  the  multivoice  discrete  undecimated  wavelet 
transform. 

[rv;j}k  =  [(Uigv)*aj]k. 

A  block  diagram  of  one  stage  of  the  transform  is  shown  in  Figure  3.3.  One  can  show 
that  (r^vs)j,jfc  and  (R^vs)jtk  can  be  computed  using  the  above  algorithm  by  choosing 
the  initializations  in  Equations  3.30  and  3.31,  respectively. 

3.4  The  Inverse  Discrete  Wavelet  Transform 
Although  the  invertibility  of  an  orthogonal  or  biorthogonal  wavelet  transform 
is  relatively  straight  forward,  the  same  cannot  be  said  with  regards  to  a  nonorthogonal 
wavelet  transform.  The  standard  inversion  procedure  for  a  nonorthogonal  wavelet 
transform  is  given  by  Equation  2.24.  However,  in  order  to  use  this  procedure  in 
practice  the  following  approximations  are  involved  [Dau90,  Dau92.  She96]:  (1)  The 
wavelet  coefficients  themselves  are  approximated,  (2)  The  duals  are  replaced  by  the 
analyzing  wavelet,  and  (3)  The  infinite  sums  are  replaced  by  partial  (finite)  sums. 
Each  of  these  approximations  can  be  justified  as  follows: 

1.  The  wavelet  coefficients  can  be  approximated  by  using  ip'(t)  as  the  analyzing 
wavelet  instead  of  %j){t)  (see  Section  3.2).  It  is  not  difficult  to  argue  that  if  tjj'{t) 
and  ip{t)  are  "close"  then  the  sequences  (s,^-)Jk)  and  (s,ipjik)  should  be  "close" 
as  well  (i.e.,  (C0s)Jifc  w  ((ys)j\fc  =  [cj{s)]k). 

2.  Replacing  the  duals  by  the  analyzing  wavelet  is  a  valid  approximation  as  long 
as  the  analyzing  wavelet  generates  a  "snug  frame"  (see  Section  2.2.1). 


3.  Replacing  the  infinite  sums  by  partial  sums  is  justified  by  the  time-frequency 
localization  properties  of  the  wavelet  transform.  It  can  be  shown  that  if  a 
signal  s(t)  is  essentially  localized  to  some  rectangle  in  the  time-frequency  plane 
[-T.  T]  x  ([-fii.  -n0]U[flo,  fii])  then  the  nonzero  wavelet  coefficients  lie  within 
or  close  to  that  rectangle  [Dau92]. 

The  above  inversion  procedure  can  be  written  as  follows 

+  00  +00 

»(*)  =  E  £  (c**)*.*^*) 

j=—oo  k=—oo 

+  O0  +O0 

j=—oo  k=-oc 

*   A   £   c,,^*)  (3-34) 

in  S 

where  5  =  {(j,  fc)  :  ro(^,*)  €  [-T, T]  and  ±  m(4*)  G  (h^i,  -«o]  U  [fi0,«i])}  and 
A'  is  a  constant  determined  by  the  dual  approximation.  In  most  practical  cases 
m{iphk)  =  2jk  and  m{^k)  =  2~Jm(^)  (see  Equations  2.14  and  2.18). 

In  practice  the  finest  scale  is  set  to  1  and  the  wavelet  coefficients  are  computed 
only  for  j  >  0  by  using  the  initialization  in  Equation  3.31.  Including  this  restriction 
in  Equation  3.34  yields 

L  1 

S(t)^EEc*w  (3-35) 

j=0  k 

where  L'  <  log2  + 1  and  the  sum  in  A;  is  implicitly  assumed  to  be  restricted  to 

k  £  [-|r>Ir]-  Although  these  approximations  are  indeed  valid  the  above  inversion 
procedure  has  several  drawbacks.  The  most  noticeable  problems  are  [She96]:  (1) 
It  does  not  provide  a  good  inverse  for  narrow-band  wavelets  (loose  frame  bounds) 
and  (2)  While  a  partial  expansion  works  quite  well  for  many  (sufficiently  oscillatory) 
signals,  it  fails  to  achieve  good  accuracy  or  requires  an  excessive  number  of  scales  for 
others.  An  alternative  solution  that  addresses  the  problems  of  the  above  inversion 
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procedure  is  presented  in  [She93].  Unlike  Equation  3.35  the  solution  in  [She93]  is 
based  on  the  undecimated  wavelet  transform  coefficients.  In  this  section  we  review 
the  main  results  of  this  work. 

Consider  L  stages  of  the  discrete  wavelet  transform  (DYYT).  Then,  the  DWT 
maps  discrete  signals  R(Z)  into  CL+1  (Z) 

DWT  (  n  r      1    „  1 

s  ->  [Wj  :  j  =  0, . . .  ,  L  -  1;  sL) 

where  (Z)  indicates  a  space  of  functions  on  the  integers  and  Wj  refers  to  the  coeffi- 
cients of  either  the  decimated  DWT,  Cj,  or  the  undecimated  DWT,  ry  For  the  above 
transform  one  has  that  [She93] 

1.  Under  very  weak  regularity  conditions,  the  norm  of  sL  goes  to  zero  as  L  goes 
to  00. 

2.  For  the  transform  to  be  nonsingular  for  finite  L,  one  must  include  the  smoothed 
signal  sL  as  well  as  the  wavelet  coefficients  Wj. 

3.  The  transform  is  linear  (although  not  time  invariant  in  the  decimated  case)  and 
may  be  represented  by  a  matrix  A. 

4.  For  /  and  g  finite,  the  transform  is  locally  finite  but  the  matrix  itself  is  infinite 
because  convolution  acts  on  arbitrarily  long  signals. 

5.  A  true  inverse  does  not  exist  because  the  image  of  the  transformation  is  a 
proper  subset  of  CL+1(Z).  However,  the  transform  is  generally  injective  so 
A]A  is  nonsingular,  and  one  may  invert  objects  in  the  range  space  CL+1  (Z)  by 
using  the  pseudo-inverse  [She93] 

P=  [A^A)~XA]. 

More  generally,  if  B  is  a  matrix  such  that  BA  is  nonsingular  then 

Q  =  (BA)~lB  (3.36) 
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Figure  3.4:  Block  diagram  of  one  stage  of  the  inverse  discrete  undecimated  wavelet 
transform. 

is  also  a  left  inverse.  However,  the  pseudo-inverse  is  unique  to  the  extent 
that  it  is  an  orthogonal  projection  of  CL+1(Z)  onto  the  image  of  the  DWT 
[She93]3.  (Compare  this  analysis  to  the  discussion  on  frames  in  Section  2.2.1.) 
Although  the  pseudo-inverse  is  a  natural  choice,  other  issues  are  involved  such 
as  computational  complexity,  or  the  the  error  criterion  not  being  a  metric  on 
CL+1(Z)  [She93]. 

An  alternative  approach  is  to  invert  a  single  stage  of  the  DWT,  thereby  in- 
verting the  entire  transform.  For  the  undecimated  DWT,  it  is  sufficient  to  find  filters 
/  and  g  such  that 

f*f  +  g*g  =  S.  (3.37) 

For  the  multivoice  case  we  have 

/*/  +  X>*0*  =  <5  (3-38) 

ti=i 

Thus  the  inverse  filter  bank  of  Figure  3.4  (Figure  3.5)  is  a  left  inverse  for  the  un- 
decimated DWT  of  Figure  3.2  (Figure  3.3)  provided  that  Equation  3.37  (Equation 
3.38)  is  satisfied.  In  fact,  as  we  will  see  shortly,  if  one  consider  a  more  general  class 
of  filters  than  those  satisfying  Equation  3.37  (Equation  3.38  in  the  multivoice  case), 
the  inverse  filter  bank  provides  a  unified  framework  from  which  to  study  various 

3Note  that  this  holds  for  non-Euclidean  metrics  on  Ci+1  (Z)  provided  that  A*  is  replaced  by  the 
adjoint  transform  relative  to  the  metric  [She93]. 
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Figure  3.5:  Block  diagram  of  one  stage  of  the  inverse  multivoice  discrete  undecimated 
wavelet  transform. 

inverses.  From  now  on,  we  will  adopt  the  terminology  DWT  to  refer  to  the  undeci- 
mated discrete  wavelet  transform  with  or  without  voices  and  inverse  discrete  wavelet 
transform  or  DWT"1  [She93]  to  refer  to  any  inverse  filter  bank  whether  or  not  the 
filters  satisfy  Equation  3.37  (Equation  3.38  in  the  multivoice  case).  Although  under 
this  framework  the  inverse  may  not  be  exact,  for  it  to  be  useful  the  filters  should  be 
chosen  to  provide  an  approximate  left  inverse  to  the  DWT.  An  exact  left  inverse  is 

achieved  by  iterating  (see  Neumann  inverse  in  Section  3.4.3). 
3.4.1    Approximations  to  the  Continuous  Inverse 

The  approach  used  in  [She93]  to  find  suitable  inverses  is  to  use  the  discretiza- 
tion of  the  continuous  inverse  in  Equation  2.21  as  a  guideline  to  provide  metrics  and 
appropriate  filter  normalizations  for  the  DWT"1.  Below  we  include  the  results  for 
the  undecimated  wavelet  transforms  with  voices.  The  decimated  wavelet  transform 
with  voices  is  treated  as  well  to  draw  some  useful  comparisons. 

From  Equation  2.21  we  have  that 


2    f+oc  f+oc  da 
s(t)  =  —  /      (Wts)(a,b)iPa;b(t)  —db. 


(3.39) 


In  the  case  of  A7  voices,  discretizing  the  integral  for  the  decimated  (a  =  2J2("~1)//V, 
b  =  Vk)  and  undecimated  (a  =  2j2^v~x^N ,  b  =  k)  case  one  obtains  [She93] 

N 


8(t) 


2  In  2 


1 

22(c^s)j,k2iv_l)/N^vj,k(t) 


(3.40) 


12 


and 

respectively;  where 

*(t)  =  v^W^  (20^)  ' 

V  £ 

and 

=  ^  H-ir)  ■ 

Notice  that  the  above  normalization  for  ipv{t)  is  not  the  same  as  the  one  in  Sec- 
tion 2.2.3.  The  normalization  used  here  follows  from  the  sampling  of  x/)a;b{t)  with 
a  =  2J2(t'_1)/yv  and  b  =  2jk  in  the  decimated  case  (or  b  =  k  in  the  undecimated 
case).  When  discretizing  Equation  3.39  one  requires  the  above  normalization  since 
the  continuous  inverse  is  defined  in  terms  of 

By  incorporating  the  square  root  of  2{v\)/N  in  the  above  normalization  oi^)v{t) 
one  can  eliminate  this  factor  from  Equations  3.40  and  3.41.  This  yields  the  normal- 
ization used  in  Section  2.2.3.  In  such  a  case,  expressions  equivalent  to  Equations  3.40 
and  3.41  are  given  by 

212^ 

v      v=l  j,k 

and 

JV 

^      v=l  j,k 

respectively.  The  normalizations  used  in  this  section  are  favored  by  [She93]  be- 
cause the  squares  of  all  four  types  of  wavelet  coefficients  (decimated/undecimated 
with/without  voices)  represent  power  per  Hz. 
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In  order  to  relate  the  above  approximations  to  the  discrete  wavelet  transform 
Shensa  [She93]  discretizes  "Plancherel"  formula  for  wavelets  (see  Equation  2.11). 
Restricting  a  >  0  we  have  that 

\\e(t)\\l  =  7r        /     l(>V)M)|2  -zdb. 
Discretizing  for  the  decimated  and  undecimated  case  with  voices  one  obtains 

^      v=l  j,k 

and 

IW«)II1«^EEm<^wI(r*.^I2-  <3-43) 

respectively.  The  goal  is  to  use  the  above  discretization  as  a  guideline  for  obtain- 
ing an  expression  equivalent  to  Plancherel's  formula  for  the  discrete  undecimated 
wavelet  transform.  Before  proceeding  in  this  direction  we  provide  an  interpretation 
for  the  additional  i  factor  appearing  in  Equations  3.41  and  3.43  when  compared  with 
Equations  3.40  and  3.42,  respectively.  As  we  will  see  shortly,  this  will  prove  useful  in 
establishing  an  expression  equivalent  to  Plancherel's  formula  for  the  discrete  wavelet 
transform. 

It  is  not  difficult  to  see  that  for  scales4  coarser  than  1  (i.e.,  j  >  0)  the  undec- 
imated wavelet  transform  has  a  denser  sampling  than  the  decimated  wavelet  trans- 
form. In  fact,  for  any  j  >  0  there  are  2j  coefficients  in  the  undecimated  wavelet 
transform  for  every  coefficient  in  the  decimated  wavelet  transform.  Similarly,  for 
scales  finer  than  1  (i.e.,  j  <  0)  we  have  the  opposite  behavior;  for  any  j  <  0  there 
are  2~J  coefficients  in  the  decimated  wavelet  transform  for  every  coefficient  in  the 
undecimated  wavelet  transform.  Therefore,  the  additional  ^  factor  in  Equations 
3.41  and  3.43  can  be  "viewed"  as  balancing  the  redundancy  at  scales  coarser  than 

4 Although  a  =  2j2{v~1]/N  is  the  scale  parameter,  in  this  discussion  we  use  "scale"  to  refer  to  the 
term  2J  . 
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1  and  compensating  for  the  conciseness  at  scales  finer  than  1.  Now  consider  the 
undecimated  DWT  of  a  discretized  signal  s  (e.g.,  see  Equation  3.31).  Assuming  the 
approximation  {Rtvs)j:k  «  [rv;j]k  we  write  the  following  "equivalent"  to  Plancherel's 
formula  in  Equation  3.43  for  the  discrete  wavelet  transform 


s 


.         N    +oo  -. 

|2_1E±VV  '  Mr-  II2 

v      u=l  j=0 


If  there  are  only  L  stages  (levels)  then  sL  "replaces"  the  coefficients  r„;j  for  j  >  L  as 
follows 


Ci,iAfZ^Z^2J2l"-1)/N"  °a  2' 


where  the  factor  ^  is  added  to  balance  the  redundancy  of  the  undecimated  wavelet 
transform  at  the  Lth  stage  (For  a  different  view  on  the  above  "derivation"  the  reader 
should  consult  [She93]).  When  the  approximation  in  Equation  3.44  holds  it  implies  a 
relationship  between  the  normalizations  of  the  filters  /  and  gv.  Consider  the  output 
of  a  one-stage  filter  bank  where  the  input  is  an  impulse  function.  Then,  Equation 
3.44  implies 

 1  IIq  l!2  +  -!l/ll2«i-  (3-45) 

^  v—\ 

Notice  that  Equations  3.44  and  3.45  are  only  meant  to  be  conceptual  guidelines  and 

will  be  used  in  the  next  section  to  define  an  appropriate  expression  for  the  energy  of 

a  (nonorthogonal)  discrete  wavelet  transform. 
3.4.2    Discrete  Inverses 

In  this  section  we  review  the  conditions  for  finite  energy  and  boundedness 
of  the  DWT  and  the  DWT"1  as  the  number  of  stages  L  becomes  infinite.  These 
properties  are  important  to  the  extent  that  they  reflect  directly  on  the  numerical 
stability  of  the  algorithms  themselves  [She93].  Although  in  practice  it  is  not  possible 
to  compute  an  infinite  number  of  stages,  the  poor  numerical  behavior  of  unbounded 
transforms  tends  to  persist  for  finite  L  [She93]. 
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Based  on  Equations  3.44  and  3.45  one  may  define  the  energy  of  a  (nonorthog- 
onal)  discrete  wavelet  transform  as  follows  [She93]:  Suppose  there  exists  cE  >  0  suc  h 
1  hat 

^  (^E  J^w)  +  ^II/H2  =  1- 

Then  define  the  energy  by  [She93] 

j=0  v-l 

This  is  to  be  interpreted  to  mean  that  {y/cErv;j,  aL)  is  an  energy  density  for  the  DWT 
and  the  total  energy  is  given  by  applying  the  metric  given  by  the  (LN+l^dimensional 
diagonal  matrix  with  elements  {2^2^v~l^N ,  2~L}  [She93]. 

It  follows  from  Equation  3.46  that  for  L  =  oo  the  discrete  wavelet  transform 
is  bounded  and  has  a  bounded  inverse  if  and  only  if  there  exist  0  <  A  <  B  <  oc  such 
that  [She93] 

A\\s\\2  <  E(r)  <  B\\s\\2. 

The  following  is  a  sufficient  condition  for  it  to  be  satisfied  in  the  case  of  a  single  voice 
[She93] 

Q  <  A       \c.E\g{en\     <B<oc  f()r  all  w<  (3.47) 
"  l-i|/(e*»)r  " 

At  u)  =  0,  Equation  3.47  is  also  a  necessary  condition  [She93].  Hence,  one  must 
have  \  <  1-  Since  /(l)  =  s/2  (see  Equation  3.27),  the  condition  in  Equation 

3.47  requires  J2k  9k  =  0  that  is'  that  9  be  a  highpass  filter  [She93].  (In  the  case  of 
multiple  voices  all  gv  must  be  highpass  filters  [She93].)  Boundedness  of  the  DWT"1 
for  arbitrary  L  also  imposes  the  constraint  [She93] 
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Similarly,  if  the  composition  of  the  inverse  and  forward  transform  is  to  be  bounded 
for  arbitrary  L  [She93] 


<  1. 


A  sufficient  condition  for  the  composition  of  DWT  1  and  DWT  to  be  a  bounded 
transform  may  also  be  derived  [She93]  and  is  given  by 

/V") I '(O  +  C~l  11(^)9(^)1  <  1  for  all  w, 


ith  1  <  C  <  oc.  Such  a  C  exists  provided 

\9(ej")9(ejw)\ 


<  1  and 


1  - 


<  C  for  all  w. 


f(ei")f(e>") 

Although  it  would  be  desirable  that  the  DWT"1  be  an  exact  left  inverse  (see 
Equations  3.37  and  3.38)  in  practice  this  leads  to  filters  that  are  too  long  for  most 
applications  [She93].  As  mentioned  before,  an  alternative  approach  is  to  provide  an 
approximate  left  inverse  to  the  DWT  and  refine  it  by  iterating  (see  Section  3.4.3). 
In  [She93],  three  filter  pairs  that  satisfy  the  above  necessary  conditions  and  provide 
useful  approximations  are  studied.  Here  we  consider  the  filter  pairs  that  are  relevant 
to  our  study:  (1)  The  double  integral  type  filters  /  =  S/y/2  and  g  =  cg\  and  (2) 
The  adjoint  filters  /  =  f /2  and  g  =  cg] .  The  filter  pair  not  considered  here  can  be 
used  only  if  ip(t)  is  analytic  or  is  real.  The  reader  can  verify  that  the  wavelet 
used  in  this  study  (see  Chapter  4)  does  not  satisfy  any  of  the  above  conditions. 

By  considering  the  output  of  the  DWT-1  in  the  case  of  a  single  voice  [She93] 

s = xin  [{u>f)  *]  (uis)  *  ri + n  [{uif)  *]  sl  (3-48) 

j=0  i=0  j=o 

one  may  interpret  the  above  filters  as  approximations  to  the  continuous  inverse  for- 
mula in  Equation  3.39.  Furthermore,  for  the  adjoint  filters  it  can  be  shown  that  the 
DWT^1  is  indeed  the  adjoint  of  the  DWT  under  the  energy  metric  in  Equation  3.46 
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[She93].  Below  we  include  the  interpretation  for  the  double  integral  type  filters.  For 
the  adjoint  filters  consult  [She93]. 

For  the  double  integral  type  filters  one  can  show  that  [She93] 

[S]o   =  ^4y[(^V)*r,]0  +  DC 
j=o  v" 

j=0  VZJ  k 

=   ^EEtcA^(-A:)  +  DC 
j=o   fc  v" 

=  ^EEf^(°)+DC  (3-49) 

j=0  fc 

where  DC  refers  to  the  last  term  of  Equation  3.48.  This  "coincides"  up  to  the  DC  term 
with  the  standard  approximation  in  Equation  3.35  at  time  t  =  0.  By  incorporating 
the  dependency  on  s  in  Equation  3.49  we  have  that 

L-l 
j=o  fc 

Since  the  undecimated  transform  and  its  inverse  are  time  invariant  one  may  write 
[She93] 

i~<  =  cEEfe^-"5)^^'0'  +  DC-  (3-50) 

3=0  k 

It  is  not  difficult  to  show  that  the  above  equation  represents  a  discrete  approximation 
to  the  double  integral  formula  in  Equation  3.39.  Consider 


2      /-oo    /»oo  ^a 

s(t  +  T)  =  (T.Ts)(t)   =   77-/     /    (W0(r_Ts))(a,6)^o;6(t)  —db 

W  J -oo  J0  ° 
2      roo    poo  ^Q 

=   77-/     /    (Wv,s)M  +  r)^ai6(f)  —  db. 

J -oo  Jo  a 


At  t  =  0  this  vields 


2  r°°  r°°  da 

s(t)  =  —        /    (Wts)(a,b  +  r)ipa,b(0)^db. 

^tp  J -oo  JO  a 
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Diseretizing  (see  the  discussion  about  the  restrictions  on  j  and  k  in  relation  to  Equa- 
tion 3.35) 

91  o1'"1 

s(n)  ~  E^K2^  2'k  +  n)^(°)- 

^     j=0  k 

But  (W*s)(2*,2»ib+n)  =  (WV)(T_ns))(2^'!  2'k)  =  (Ct(T_ns))jM  w  [Cj(T_„«)]*.  There- 
fore, 

*(«)  *     E  E^T-»s)^(°)-  (3-51) 

This  "coincides"  up  to  the  DC  term  with  Equation  3.50.  Notice  that  Equation  3.50 
should  be  a  better  approximation  since  the  DC  term  replaces  the  energy  lost  in 
truncating  the  expansion  at  j  =  L  -  1. 

The  above  formulation  is  extended  to  jV  voices  by  considering  the  following 
filters  [She93]:  (1)  The  double  integral  type  filters 

/  =  S/V2  and  gv  =  ^fffOl  (3-52) 

and  (2)  The  adjoint  filters 

f  =  f^2  aixd  gv  =  y^ol  (3-53) 

For  additional  details  the  reader  should  consult  [She93]. 

Finally,  we  briefly  review  the  selection  of  the  parameter  c  for  the  above  filters. 
Although  there  are  several  approaches  proposed  in  [She93]  we  only  review  the  most 
successful  ones  here:  (1)  Setting  the  energy  of  the  impulse  response  of  a  single  stage 
of  the  composition  of  DWT-1  and  DWT  to  be  equal  to  one,  i.e., 

\\ct,^yNsU9v  +  f*f\\2  =  ^  (3-54) 


17=1 

and  (2)  Setting 

21n2 


<■  = 


(3.55) 

cviv 

by  analogy  to  the  continuous  case.  (Compare  Equations  3.50  and  3.51.  The  normal- 
ization 1/N  comes  from  the  generalization  to  multiple  voices.) 
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3.4.3    Neumann  Inverse 

Given  a  discrete  wavelet  transform  A  and  an  approximate  inverse  C,  an  exact 
left  inverse  (see  Equation  3.36)  can  be  computed  by  using  the  Neumann  inverse  for 
linear  operators.  Here,  A  and  C  are  the  transforms  DWT  and  DWT"1,  respectively. 
Let  S  be  a  linear  operator  in  a  Hilbert  space  then  the  Neumann  inverse  is  given  by 

oo 

i,  (i 

where  ||X  -  fiS\\  <  1  to  insure  convergence  [Dau92,  She93].  (Notice  this  procedure 
is  also  used  to  compute  the  duals  in  Equation  2.32.)  Equation  3.56  implies  that  the 
inverse  of  CA  is  given  by 

00 

k—0 

Therefore,  given  a  signal  s  we  have  that 

oo 

s  =  ^YJ{I-^CA)kCAs. 

Truncating  after  N  +  1  terms  we  have  that 

N 

sN  =  -  fiCA)kCAs. 

k=0 

which  can  be  computed  iteratively  as  follows 
Algorithm  3.4.1 

k  <_  o,  ek  <r-  xk  <r-  DWT"1  oDWTs  {Decomposition  followed  by  reconstruction.} 
while  k  <  N  do 
k  i-  k  +  1 

ek  ^  ek-i  _  /;  DWT"1  o  DWTe*"1 
xk  <r-  xk~l  +  ek 
end  while 

sk  ^—  fixk  {kth  left  inverse.} 
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As  pointed  out  in  [She93],  Algorithm  3.4.1  suffers  from  a  computational  draw- 
back as  there  is  a  huge  increase  in  the  support  of  the  signal  under  the  mapping 
DWT"1  o  DWT.  In  [She93],  a  heuristic  solution  that  consists  in  clipping  the  support 
of  the  iterations  to  that  of  the  DWT  of  s  is  proposed.  As  we  will  see  in  Chapter  5 
we  avoid  this  problem  by  restricting  s  to  periodically  extended  signals. 


CHAPTER  4 
THE  1-D  SINE-GABOR  WAVELET 

In  this  chapter  we  introduce  a  new  wavelet:  the  sine-Gabor  function.  We 
show  that  this  function  not  only  satisfies  the  wavelet  admissibility  condition  but 
under  certain  conditions  achieves  nearly  optimum  time-frequency  localization  and 
generates  a  frame  of  wavelets.  Furthermore,  we  show  that  there  is  a  trade-off  between 
the  tightness  of  the  frame  and  the  time-frequency  localization  properties  of  the  sine- 
Gabor  wavelet.  We  demonstrate  how  this  trade  off  can  be  overcome  by  introducing 
voices. 

4.1    The  Sine-Gabor  Wavelet 

It  is  well  known  that  the  modulated  Gaussian  or  Gabor  function  is  the  only 
function  that  achieves  optimum  time- frequency  localization.  The  question  arises  as 
to  whether  the  Gabor  function, 

Vv.(f)  =  A'Ge-t2/2"oe-^0', 

can  be  used  in  the  context  of  wavelet  analysis  and  synthesis.  Unfortunately,  this  is 
not  possible  because  the  Gabor  function  is  not  an  admissible  wavelet  (for  details  on 
the  admissibility  condition  for  wavelets  see  Section  2.2.)  Therefore,  optimum  time- 
frequency  localization  via  the  wavelet  transform  is  not  possible.  However,  a  modified 
version  of  the  Gabor  function,  know  as  the  Morlet  wavelet  [Dau92], 

1  /2 

satisfies  the  wavelet  admissibility  condition,  and  for  cr0  =  T  =  7r(2/ln2)  , 
achieves  nearly  optimum  time-frequency  localization  and  generates  a  frame  of 
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wavelets  (for  more  details  about  wavelet  frames  see  Section  2.2  and  2.2.1.)  In  fact, 
the  error  incurred  in  approximating  the  Morlet  wavelet  with  a  Gabor  function  under 
these  conditions  is  negligible  [Dau92].  Although  Gabor  functions  are  not  wavelets  in 
the  strict  sense,  the  term  "Gabor  wavelet"  is  used  in  the  literature  (e.g.,  [Bov92]).  We 
will  refer  to  Gabor  functions  as  Gabor  wavelets  provided  that  e~a^l'2  [s  negligible. 

Having  a  frame  of  wavelets,  ipj,k(t),  j,k  <E  Z,  is  a  desirable  property  because 
it  allows  to  completely  characterize  and  reconstruct  a  signal  s(t)  from  its  discrete 
wavelet  coefficients  {s(t),ipjtk{t))  (see  Section  2.2.)  In  practice  it  is  desirable  to  have 
a  frame  of  wavelets  for  which  both  analysis  and  synthesis  can  be  computed  efficiently. 

We  show  that  the  sine-Gabor  function  (imaginary  part  of  the  Gabor  function) 
achieves  nearly  optimum  time-frequency  localization  and  at  the  same  time  generates 
a  frame  of  wavelets.  Furthermore,  we  show  that  both  analysis  and  synthesis  using 
the  sine-Gabor  wavelet  can  be  computed  efficiently. 

The  fact  that  the  linear  combination  of  two  wavelets  is  also  a  wavelet  gives  us 
an  easy  way  to  show  that  the  sine-Gabor  function  is  a  wavelet.  By  taking  the  real 
and  imaginary  parts  of  the  Morlet  wavelet  we  have  that 

RepM(t)  =  KAIe-x'2^o  (cos(uV)  -  e~ao^2)  . 

and 

Im^Af(t)  =  -KMe-*2/2cr2°  sinM). 

It  follows  that  the  sine-Gabor  function  is  indeed  a  wavelet.  It  also  follows  that  the 
cosine-Gabor  function  (real  part  of  the  Gabor  function)  is  not  a  wavelet  for  otherwise 
the  Gabor  function  would  be  a  wavelet. 

The  question  as  to  whether  or  not  the  sine-Gabor  wavelet  generates  a  frame 
of  wavelets  and  the  efficient  computation  of  both  analysis  and  synthesis  using  the 
sine-Gabor  wavelet  are  addressed  below. 


53 


4.2    The  Sine  Gabor  Wavelet  Frame 

We  showed  in  Section  4.1  that  sine-Gabor  functions  are  wavelets.  In  this 
section  we  study  the  time-frequency  localization  characteristics  of  the  sine-Gabor 
wavelets  and  also  establish  that  under  certain  conditions  they  generate  wavelet 
frames. 

The  sine-Gabor  wavelet  is  given  by 

4>{t)  =  K^^hmiuot),  (4.1) 

where 

\ff0(l-c"'oWb)/ 

was  chosen  so  that  ||i/>||2  =  1-  The  time  and  frequency  resolution  of  the  sine-Gabor 
wavelet  are  given  by  the  root  mean  square  extent  of  the  function  and  its  Fourier 
transform,  respectively  (see  Equations  2.1  in  Section  2.1.)  It  follows  from  Equations 
4.1  and  2.1  that  the  time-resolution  of  the  sine-Gabor  wavelet  is  given  by 

aW=(^£^^mY\  ,4.2) 

v  '  y     2(1-6-^5)  j 

The  Fourier  transform  of  the  sine-Gabor  wavelet  is  given  by 

e-ffg(w-wo)a/2  _  e-cl(uJ+LJo)2 /2 

$M  =  Kj  -  —  • 

where 


Computing  the  root  mean  square  extent  of  t/>(u>)  as  defined  in  Equation  2.1  would  lead 
to  misleading  results.  It  can  be  verified  that  if  a  function  tp{t)  satisfies  the  wavelet 
admissibility  condition  and  is  also  a  window  function  then  its  Fourier  transform 
iJ)(uj)  must  vanish  at  the  origin  (see  Equation  2.16  in  Section  2.2.)   In  addition  if 
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10  ,0 


Figure  4.1:  Square  of  the  product  of  the  time  and  frequency  resolution  for  the  sine- 
Gabor  wavelet  as  a  function  of  a0  and  u0. 


the  wavelet  tp(t)  is  real  we  have  that  =  i>(-u),  so  that  \ip{u)\  is  a  symmetric 
function  around  the  origin.  It  is  therefore  more  suitable  to  restrict  our  attention  to 
positive  frequencies  for  this  class  of  wavelets  and  use 


(o;-m(^))2|^(u;)|2  dw 


for  the  frequency-resolution,  where 

ra(V0  = 


-f 


uo\x[}(u})\2dLO. 


It  is  easy  to  verify  that  the  sine-Gabor  wavelet  satisfies  all  the  conditions  mentioned 
above  and  that  its  frequency-resolution  is  given  by 

1  +  2<jIlo'1  -  e-ffowo 


1/2 


where 


2(T2(1  -  e-^o) 

uQ  evi{aou0) 


(4.3) 


(4.4) 


Figure  4.1  shows  the  square  of  the  product  of  the  time  and  frequency  resolution  for 
the  sine-Gabor  wavelet  as  a  function  of  a0  and  uj0.  For  a0uJo  >  5/2,  the  sine-Gabor 
wavelet  is  within  1%  of  the  optimum  time-frequency  localization  lower  bound  given 
by  the  uncertainty  principle,  a2{il>)  cr2(^)  >  1/4  (for  more  details  see  Equation  2.8  in 
Section  2.1.) 
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Next,  we  establish  that  the  sine-Gabor  function  does  indeed  generate  a  frame 
of  wavelets.  General  conditions  on  a  wavelet  ip(t)  under  which  a  frame  is  obtained 
and  estimates  for  the  frame  bounds  are  discussed  in  Section  2.2.2.  A  wavelet  tp(t) 
generates  a  frame  of  wavelets  if  the  following  condition  is  satisfied 

<  C|u;|a(l  +  M)-0  (4.5) 

with  a  >  0,  /3  >  a  +  1. 

The  magnitude  of  the  Fourier  transform  of  the  sine-Gabor  wavelet  may  be 
written  as  follows 

<  ^(2a02u;ok|)1/2e-ffo(M-"o)2/2 

(2aga;ol^l)1/2 
-     2   l  +  ^da;!-^)2^' 

where  the  first  and  second  inequalities  are  established  by  using  1  -  e'^  <  u;7,  V^'  > 

0.  0  <  7  <  1,  and  e~Jil2  <  1+j,/2,  Vcj,  respectively.  It  follows  that  the  sine-Gabor 

wavelet  satisfies  Equation  4.5  with  a  =  0.5  and  0  =  1.85  for  some  C  sufficiently  large. 
4.2.1    Frame  Bounds 

Table  4.1  and  Figure  4.2  show  the  estimated  frame  bounds  and  sine-Gabor 
wavelet  for  various  values  of  a0  and  Wo-  The  results  in  Table  4.1(a)  show  that  it 
is  possible  to  obtain  a  frame  of  wavelets  by  using  the  sine-Gabor  wavelet  and  at  the 
same  time  achieve  nearly  optimum  time-frequency  localization,  a2(ip)a2(ip)  =  0.2525. 
Table  4.1(b)  shows  that  by  trading  off  time-frequency  localization,  a2(tp)a2{^)  = 
0.3297,  the  frame  bounds  can  be  made  tighter.  As  mentioned  in  Section  2.2,  snug 
frame  bounds  are  a  desirable  property  because  the  dual  wavelet  can  be  approximated 
by  jtb^(<)-  The  error  between  the  original  signal  and  its  reconstruction  using  the 
above  approximation  is  given  by  ^| I/H2,  wnere  r  =  B/A  -  1  (see  Section  2.2.1  for 
more  details.)  For  the  sine-Gabor  wavelet  with  aQ  =  1,  u>0  =  1,  the  error  is  2.55%  for 
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Table  4.1:  Frame  bounds  for  the  sine-Gabor  wavelet  with  (a)  a0  =  5/2,  up  =  1, 
(b)  a0  =  1,  u)q  =  1.  and  (c)  a0  =  1.0657,  u0  =  0.0299.  The  elementary  dilation 
parameter  q0  =  2  in  all  eases. 

(a)  Frame  bounds 
with  a0  =  5/2,  ujq  =  1. 
=  bp        A  B  BjT 

1.00  4.3873  5.6329  1.2839 
1.25  3.5099  4.5063  1.2839 
1.50  2.9216  3.7586  1.2865 
1.75  2.4089  3.3169  1.3769 
2.00  1.6016  3.4085  2.1281 
2.25    0.3963    4.0572  10.2385 


(b)  Frame  bounds 
with  cto  =  1,  u/'o  =  1. 
~    ~A       ~B~  B/A 
0.50    8.5076    8.6988  1.0225 
0.75    5.6715    5.7995  1.0226 
1.00    4.1918    4.4114  1.0524 
1.25    2.8464    4.0362  1.4180 
1.50    1.2635    4.4720  3.5395 


(c)  Frame  bounds  with 
a0  =  1.0657,  u0  =  0.0299. 
b0        A  B  BjA 

0.75  7.2030  7.3285  1.0174 
1.00  5.3999  5.4988  1.0183 
1.25  4.2456  4.4733  1.0536 
1.50  3.1604  4.1054  1.2990 
1.75  1.9203  4.3075  2.2432 
2.00    0.6524    4.7969  7.3523 


Table  4.2:  Frame  bounds  for  the  first  derivative  of  a  Gaussian  with  a0  =  1.0657. 


bp        A  B  B/A 

0.75  7.2043  7.3298  1.0174 

1.00  5.4008  5.4997  1.0183 

1.25  4.2465  4.4739  1.0536 

1.50  3.1615  4.1055  1.2986 

1.75  1.9218  4.3071  2.2412 

2.00  0.6541  4.7962  7.3328 
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Figure  4.3:  First  derivative  of  a  Gaussian  (left)  and  the  magnitude  of  its  Fourier 
transform  (right)  with  a0  =  1.0657. 
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o0  =  2, 60  =  1-  The  parameters  for  the  sine-Gabor  wavelet  used  in  Table  4.1(c)  were 

obtained  by  using  a  Simplex  search  method  to  find  the  minimum  of  the  difference 

between  B  and  .4  as  a  function  of  a0  and  u0.  The  error  in  this  case  is  0.90%  for 

«o  =  2,  b0  =  1.  This  reduced  error  is  at  the  expense  of  time-frequency  localization. 

a2(ip)a2(ip)  =  0.3401.  An  important  observation  follows  from  this  last  example:  For 

cr0  =  1.0657,  wo  =  0.0299  the  sine-Gabor  wavelet  is  nearly  equal1  (up  to  a  sign)  to 

the  first  derivative  of  a  Gaussian  with  cr0  =  1-0657.  Therefore,  the  first  derivative  of 

a  Gaussian  (Canny's  approximation  to  the  optimal  step  edge  detector  [Can86])  can 

be  used  to  completely  characterize  and  reconstruct  a  signal  from  its  discrete  wavelet 

coefficients.  Figure  4.3  shows  the  first  derivative  of  a  Gaussian  and  the  magnitude 

of  its  Fourier  transform.  Table  4.2  shows  the  frame  bounds  for  the  first  derivative 

of  a  Gaussian  with  a0  =  1.0657.  The  advantage  of  the  sine-Gabor  function  over  the 

first  derivative  of  a  Gaussian  is  that  by  changing  the  product  of  a0iv0,  one  is  able  to 

trade-off  between  Canny's  criteria  for  designing  optimal  step  edge  detectors,  that  is, 

the  product  of  signal-to-noise  ratio  and  localization,  EA,  and  the  multiple  response 

constraint,  r  [Meh92].  This  is  a  property  of  Canny's  optimal  step  edge  detector  that 

is  not  possible  to  achieve  with  the  first  derivative  of  a  Gaussian.  Table  4.3  shows  the 

frame  bounds  for  the  sine-Gabor  wavelet  for  0.75  <  a0u0  <  3.75.    Figure  4.4  shows 

the  frequency  coverage  of  Equation  2.34  for  the  sine-Gabor  wavelet.  Notice  that  this 

equation  is  in  fact  the  dominant  term  in  the  computation  of  the  frame  bounds  in  this 

case  (compare  Figure  4.4  and  Table  4.1)  and  that  by  trading  off  the  bandwidth  of 

the  wavelet  one  is  able  to  achieve  a  tighter  frame. 
4.2.2    Multi voice  Frame  Bounds 

In  the  previous  section  we  showed  that  the  sine-Gabor  wavelet  does  indeed 
generate  a  wavelet  frame  but  that  it  is  not  possible  to  generate  a  snug  sine-Gabor 
!It  is  assumed  that  both  functions  have  the  same  L2-norm. 


59 


6 


5 


4 


■  »  /  \ 


10  12  14  16 


(b) 


Figure  4.4:  Frequency  coverage  of  ^  r0(2Ju;) 


for  the  sine-Gabor  wavelet  over 


j  =  -oo 


four  octaves  (solid)  and  significant  terms  |v>(2Ju;)  (dashed)  contributing  to  the  sum 
over  the  range  shown,  (a)  a0  =  5/2,  u0  =  1,  (b)  a0  =  1,  u)q  =  1. 
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Table  4.3:  Frame  bounds  for  the  sine-Gabor  wavelet  for  0.75  <  aQuQ  <  3.75.  The 
elementary  dilation  parameter  a0  =  2  and  the  elementary  translation  step  b0  =  1  in 
all  cases. 


A 


B 


B/A 


3 

75 

2.5933 

6.8524 

2.6423 

3 

50 

2.9659 

6.5120 

2.1956 

3 

25 

3.3385 

6.2107 

1.8603 

3 

00 

3.7031 

5.9576 

1.6088 

2 

75 

4.0546 

5.7629 

1.4213 

2 

50 

4.3873 

5.6329 

1.2839 

2 

25 

4.6920 

5.5627 

1.1856 

2 

00 

4.9418 

5.5248 

1.1180 

1 

75 

5.0866 

5.4618 

1.0738 

1 

50 

5.0580 

5.2931 

1.0465 

1 

25 

4.7875 

4.9433 

1.0325 

1 

00 

4.1918 

4.4114 

1.0524 

0 

75 

2.8851 

4.0837 

1.4154 

10  12  14  16 


N  +oo 


for  the  multivoice  sine-Gabor 


Figure  4.5:  Frequency  coverage  of  ^  ^  j^^u;) 

t)=l  j=— 00 

2 

wavelet  over  four  octaves  (solid)  and  significant  terms  ■0„(2Ju;)  (colored)  contribut- 
ing to  the  sum  over  the  range  shown  for  N  =  4,  a0  =  5/2,  uq  =  1.  Colors  are  used 
to  denote  the  four  voices  within  each  octave. 
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Table  4.4:  Frame  bounds  for  the  sine-Gabor  wavelet  with  a0  -  5/2,  ui0  =  1,  (a)  N  - 
2,  (b)  N  =  3,  and  (c)  N  =  4.  The  elementary  dilation  parameter  a0  =  2  in  all  cases. 


(a)  Frame  bounds  with  N  —  2. 


b0 

A 

B 

B/A 

0.25 

40.0790 

40.1998 

1.0030 

0.50 

20.0395 

20.0999 

1.0030 

0.75 

13.3597 

13.3999 

1.0030 

1.00 

10.0197 

10.0500 

1.0030 

1.25 

8.0158 

8.0400 

1.0030 

1.50 

6.6765 

6.7033 

1.0040 

1.75 

5.6274 

5.8410 

1.0379 

2.00 

4.4175 

5.6174 

1.2716 

2.25 

2.8891 

6.0308 

2.0874 

(b)  Frame  bounds  with  N  -  3. 

bQ 

.4 

B 

B/A 

0.25 

60.2085 

60.2097 

1.0000 

0.50 

30.1042 

30.1049 

1.0000 

0.75 

20.0695 

20.0699 

1.0000 

1.00 

15.0521 

15.0524 

1.0000 

1.25 

12.0417 

12.0420 

1.0000 

1.50 

10.0314 

10.0383 

1.0007 

1.75 

8.5028 

8.6998 

1.0232 

2.00 

6.9229 

8.1294 

1.1743 

2.25 

5.0273 

8.3525 

1.6614 

(c)  Frame  bounds  with  N  =  4. 

h 

A 

B 

B/A 

0.25 

80.2788 

80.2788 

1.0000 

0.50 

40.1394 

40.1394 

1.0000 

0.75 

26.7596 

26.7596 

1.0000 

1.00 

20.0697 

20.0697 

1.0000 

1.25 

16.0558 

16.0558 

1.0000 

1.50 

13.3765 

13.3831 

1.0005 

1.75 

11.3686 

11.5682 

1.0175 

2.00 

9.4022 

10.6675 

1.1346 

2.25 

7.0938 

10.7459 

1.5148 
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wavelet,  frame  and  at  the  same  time  obtain  nearly  optimum  time-frequency  local- 
ization. This  limitation  can  be  overcome  by  the  introduction  of  voices  (see  Section 
2.2.3).  As  mentioned  before  voices  not  only  aid  in  obtaining  a  tighter  frame  but  do 
so  without  giving  up  too  much  freedom  in  choosing  the  wavelet  or  its  bandwidth.  In 
addition,  voices  or  suboctave  sampling  is  a  desirable  property  for  applications  in  the 
area  of  signal  analysis  and  time-frequency  localization  [She96].  Table  4.4  shows  the 
frame  bounds  for  the  multivoice  sine-Gabor  wavelet  using  the  estimates  described  in 
Section  2.2.3  and  fractionally  dilated  versions  of  a  single  wavelet  as  described  in  Equa- 
tion 2.35.  Notice  that  with  two  voices  {N  =  2)  the  error  between  the  original  signal 
and  its  reconstruction  using  -£gip(t)  as  the  dual  wavelet  is  0.15%  for  a0  =  2,  b0  -  1 
(see  Section  4.2.1  for  a  discussion  on  the  error.)  This  is  achieved  while  preserving 
nearly  optimum  time-frequency  localization  <r2(V>)  cr2(-0)  =  0.2525  (within  1%  of  the 
optimum).  For  N  =  3  and  N  =  4  the  error  is  negligible  for  a0  =  2,  bQ  =  1-  In  con- 
trast, the  error  with  TV  =  1  is  12.43%  for  a0  =  2,  bQ  =  1  and  the  same  time-frequency 
resolution  (see  Table  4.1(a)).  Figure  4.5  shows  the  frequency  cover  of  Equation  2.36 
for  the  sine-Gabor  wavelet  with  N  —  4. 


CHAPTER  5 

THE  I  D  DISCRETE  SINE-GABOR  WAVELET  TRANSFORM 

In  this  chapter  we  show  how  to  use  and  extend  the  approach  introduced  in 
Chapter  3  to  compute  the  forward  and  inverse  discrete  wavelet  transform  using  the 
sine-Gabor  function  as  the  analyzing  wavelet.  We  also  evaluate  the  quality  of  the 
approximating  wavelet  ip' (t)  for  the  sine-Gabor  wavelet. 

5.1    A  Trous  Filters  and  Wavelets 

Before  turning  into  the  computation  of  the  discrete  wavelet  transform  using 
the  sine-Gabor  function  as  the  analyzing  wavelet,  we  evaluate  the  quality  of  the 
resulting  approximation  tp' (t)  for  the  sine-Gabor  wavelet  given  a  suitable  a  trous 
filter  (see  Equations  3.8,  3.29  and  4.1).  Given  an  a  trous  filter,  the  Fourier  transform 
of  the  corresponding  scaling  function  is  given  by  Equation  3.25.  It  is  assumed  that 
suitable  regularity  conditions  are  imposed  so  that  the  inverse  Fourier  transform  of 
the  product  in  Equation  3.25  converges  to  a  reasonably  behaved  function  (for  more 
details  see  [Dau88,  Mal89b,  She92]).  The  wavelet  ip' (t)  is  then  given  by  combining 
Equations  3.26  and  3.29,  where  gk  are  the  samples  of  the  sine-Gabor  wavelet  tp(t)  in 
Equation  4.1,  i.e.,  gk  =  il'(-k).  Unfortunately,  it  is  not  possible  to  compute  either  the 
scaling  function  or  the  wavelet  ip'  (t)  as  this  requires  an  infinite  number  of  iterations. 
In  order  to  be  able  to  evaluate  the  quality  of  the  approximating  wavelet  we  instead 
compute  a  finite  number  of  iterations  and  compare  the  samples  of  the  approximating 
wavelet  to  the  exact  samples  of  the  sine-Gabor  wavelet.  Based  on  Equation  3.26  we 
define  the  jth  iteration  of  the  scaling  function  <j>{t)  as 

<t>u)(t)  =  Yy(DF)3]^x(Vt  -  k).  (5.1) 
k 
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It  follows  that 

<f>(t)  =  lim  4>U)  (t). 

Based  on  Equations  3.29  and  5.1  and  with  the  aid  of  Equation  3.6  we  define  the  j'th 
iteration  of  the  approximating  wavelet  if/ (t)  as 

4\3)(t)  =  YyG(DF)%^x(23t  -  k).  (5.2) 

A: 

It  follows  that 

t{>'(t)=  lim  ip'U)(t). 
From  Equations  5.1  and  5.2  it  follows  that 

fa  (£)  =  V2J[(DF)%n  (5-3) 

and 

4)  (£)  =  ^[G(DF)%,n,  (5.4) 
respectively.  One  can  show  that  the  terms  in  Equations  5.3  and  5.4  are  given  by 

[(DFy)0,n  =  [fUuf*---*uj-lf}n 

and 

[G(DF)%n  =  [f  *  Uf  *  ■  •  •  *  U^f  *  C/V]n 

respectively. 

Before  computing  the  scaling  function  and  the  approximating  wavelet  we  still 
need  to  resolve  three  issues:  (1)  Which  a  trous  filter  fk  are  we  going  to  use?,  (2)  At 
what  rate  are  we  going  to  sample  ip(t)  to  obtain  gk?,  and  (3)  How  many  samples  of 
xp(t)  are  going  to  be  used  to  define  gk? 
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Table  5.1:  Filter  coefficients  of  Lagrange  a  trous  filters  of  various  lengths  P  {2Q  -  1 
is  the  degree  of  the  corresponding  Lagrange  polynomial,  P  =  AQ  -  1). 


/ 

P  =  Z,Q  =  1 

P  =  7,  Q  =  2 

P  =  ll,  Q  =  3 

p=  15,  Q  =  4 

/o 

n/2/2 

s/2/2 

0.7071 

0.7071 

f±l 

v/2/4 

9v^/32 

0.4143 

0.4230 

f±2 

0 

0 

0 

f±3 

-V2/32 

-0.0691 

-0.0846 

/±4 

0 

0 

./  :  "» 

0.0083 

0.0169 

/±6 

0 

f±7 

-0.0017 

The  only  restriction  in  choosing  the  a  trous  filter  is  that  it  must  lead  to 
a  reasonable  behaved  scaling  function.  We  arbitrarily  adopt  the  commonly  used 
Lagrange  a  trous  filters  [She92].  Table  5.1  shows  the  filter  coefficients  of  Lagrange 
a  trous  filters  of  various  lengths.  With  regards  to  the  sampling  of  ip(t)  we  will 
use  Nyquist  criterion.  Since  the  sine-Gabor  wavelet  is  neither  bandlimited  nor  time 
limited  we  need  to  develop  criteria  for  determining  the  sampling  frequency  as  well  as 
the  number  of  samples  to  be  used.  We  use  the  mean  and  standard  deviation  of  the 
wavelet  and  its  Fourier  transform  to  establish  such  criteria.  In  the  case  of  sine-Gabor 
wavelet  we  refer  to:  (1)  a(xp)  (see  Equation  4.2)  as  the  "support1'  of  the  wavelet, 
(2)  m(0)  (see  Equation  4.4)  as  the  center  frequency  of  the  wavelet,  and  (3)  o(ip) 
(see  Equation  4.3)  as  the  "bandwidth"  of  the  wavelet.  It  is  easy  to  verify  that  for 
(j0uj0  >  5/2  these  quantities  can  be  approximated  by 

oW  -  ^, 

m(ip)  «  lu0, 

and 
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We  determine  the  maximum  frequency  for  the  sine-Gabor  wavelet  by 

Umax  =  m(i>)  +  duo{i)), 

where  the  value  du  is  an  integer  such  that  the  energy  of  in  the  interval  [-ow, 
Umax)  is  larger  than  a  certain  percentage  of  its  total  energy.  Similarly,  we  truncate 
the  sine-Gabor  wavelet  to  the  time  interval 

[-dt(r{il>),dta(i))], 

where  the  value  dt  is  an  integer  such  that  the  energy  of  ip(t)  in  the  above  interval  is 
larger  than  a  certain  percentage  of  its  total  energy. 

In  order  to  demonstrate  this  procedure  consider  computing  the  approximating 
wavelet  ip'(t)  for  the  sine-Gabor  wavelet  ip(t)  with  a0  =  5/2  and  wo  =  1.  Recall  that 
these  values  lead  to  a  wavelet  that  is  within  1%  of  the  optimum  time-frequency 
resolution.  First  we  must  choose  a  sampling  rate.  Although  the  Fourier  transform  of 
the  sine-Gabor  wavelet  is  not  bandlimited  one  may  verify  that  for  dw  =  3  more  than 
99.85%  of  the  energy  of  is  confined  to  the  interval  u  €  [-1.85, 1.85].  Therefore, 
a  sampling  frequency  u)s  >  3.7  should  suffice  in  practice.  This  implies  that  the 
sampling  period  Ts  <  1.69.  For  simplicity  we  choose  Ts  =  1.  Figure  5.1  shows  the 
sine-Gabor  wavelet  and  its  sampled  version  for  Ts  =  1.  Next  we  must  truncate  the 
number  of  samples  to  be  used.  One  may  verify  that  for  dt  =  3  more  than  99.84% 
of  the  energy  of  ip{t)  is  confined  to  the  interval  t  €  [-5.37,5.37].  Therefore,  we 
define  gk  =  ip{-k).  -6  <  k  <  6,  k  e  Z.  Table  5.2  shows  the  filter  coefficients 
gk  =  ip(-k),  -6  <  k  <  6.  k  6  Z  for  a0  =  5/2  and  ljq  =  1.  As  an  example  we 
compute  the  iterated  scaling  function  and  approximating  wavelet  with  gk  as  defined 
above  using  the  Lagrange  a  trous  filter  of  length  P  =  7.  Figures  5.2  and  5.3  show  the 
scaling  function  and  the  approximating  wavelet  for  six  iterations.  The  error  between 
^0)  i%)  and  V>  (57)  f°r  tne  tentn  iteration  is  shown  in  Figure  5.4.  Notice  that  indeed 
ip'(k)  -  i)(k)  =  0  for  -6  <  k  <  6,  k  G  Z  as  expected  (see  Equation  3.29).  Numerical 
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Figure  5.1:  Sine-Gabor  wavelet  for  a0  =  5/2  and  u0  =  1  and  corresponding  samples 
for  Ts  =  1. 


Table  5.2:  Filter  coefficients  gk  =  i'(-k)  for  cr0  =  5/2  and  u0  =  1. 
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Figure  5.2:  Iterations  1  through  3  for  (a)  The  scaling  function  0(f),  and  (b)  The 
approximating  wavelet  ip'(t).  Samples  are  shown  with  a  '+'  sign  scale  permitting. 
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Figure  5.3:  Iterations  4  through  6  for  (a)  The  scaling  function  <j){t),  and  (b)  The 
approximating  wavelet  ip  (t). 
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Figure  5.5:  Block  diagram  of  a  2-level  forward  multivoice  discrete  undecimated 
wavelet  transform. 

results  show  that  the  relative  integrated  root  mean  square  error  between  ip[10){t)  and 
ip(t),  i.e.,     '  °^(t)[|a  ,  is  0.0967c 

5.2    Forward  Discrete  Sine  Gabor  Wavelet  Transform 

In  this  section  we  address  the  implementation  issues  regarding  the  forward 
discrete  sine-Gabor  wavelet  transform.  In  particular  we  are  interested  in  computing  a 
multivoice  discrete  undecimated  wavelet  transform.  An  example  of  a  2-level  transform 
is  shown  in  Figure  5.5.  Since  our  implementation  will  be  carried  out  on  a  digital 
computer  we  need  to  constrain  input  signals  sk  to  have  finite  length.  Without  loss  of 
generality  we  will  assume  that  sk  is  defined  for  k  G  [0,  K  -  1].  Wre  take  advantage  of 
this  limitation  to  constrain  the  length  of  the  output  signals  of  the  filter  bank  to  have 
the  same  length  as  the  input.  This  is  accomplished  by  using  as  input  the  periodic 
extension  of  sk.  As  a  result  each  of  the  output  signals  of  the  filter  bank  is  periodic 
with  the  same  period  of  the  input  and  therefore  can  be  represented  with  only  K 
samples.  However,  this  limits  the  length  of  the  original  input  signal  to  be  at  least 
as  long  as  the  longest  filter  in  the  filter  bank.  The  reason  for  this  constrain  is  that 
we  want  to  preserve  the  shape  of  the  filter  kernels.  If  the  original  input  signal  were 
to  be  shorter  than  a  certain  filter  kernel  then  the  periodic  extension  of  the  signal 
would  amount  to  a  wrap  around  of  the  filter  kernel.  In  such  a  case  the  output  of 
the  filter  bank  will  bear  no  relation  to  the  wavelet  "series"  coefficients  we  are  seeking 
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to  compute  due  to  the  change  in  the  filter  kernels.  This  restriction  can  be  viewed 
as  a  limit  to  the  number  of  levels  that  can  be  computed  in  practice.  Let  M  be  the 
maximum  filter  length  amongst  the  filters  f,g},...  ,gN,  and  let  K  be  the  length  of 
the  original  input  signal,  then  the  maximum  number  of  levels  L  that  can  be  computed 
is  given  by 


It  can  be  shown  that  given  a  signal  sk,  k  £  [0,  K-l]  the  above  implementation  would 
amount  to  limiting  the  number  of  levels  to  L  or  less,  using  sk  as  input,  and  performing 
circular  convolution  instead  of  convolution  in  the  filter  bank.  One  drawback  of  this 
method  is  that  it  may  introduce  artifacts  in  the  outputs  of  the  filter  bank  as  extending 
the  input  signal  periodically  can  result  in  an  abrupt  transition  not  present  in  the 
original  signal.  This  problem  can  be  alleviated  by  using  as  input  signal  the  mirror 
extension  of  the  original  input  signal.  Given  sk,  k  e  [0,  K  -  1]  its  mirror  extension 
mk  is  given  by 


Although  from  now  on  we  will  assume  that  the  input  to  the  filter  bank  is  the  mirror 
extension  of  the  original  input  sk  we  still  limit  the  number  of  levels  to  be  computed 
to  L  in  Equation  5.5.  Before  giving  a  justification  for  maintaining  this  constrain  we 
need  to  study  the  signals  at  the  outputs  of  the  filter  bank. 

Following  a  scheme  similar  to  the  one  presented  in  [Kor96,  Opp89],  we  define 
the  following  generalized  symmetric/antisymmetric  signal  types: 

•  A  signal  sk,  k  €  [ss,  ss  +  K  -  1]  satisfying  sSs+k  =  sSt+K-i-k,  k  €  [0,  [^J]  is 
said  to  be  of  Type  I  if  K  is  odd  and  of  Type  II  if  K  is  even. 

•  A  signal  sk,  k  £  [ss,ss  +  K  -  1]  satisfying  sSs+k  =  -sSt+K-i-k,  k  E  [0,  [^J] 
is  said  to  be  of  Type  III  if  K  is  odd  and  of  Type  IV  if  K  is  even. 


log2 


K  -  1 
M  -  1 


+  1. 


(5.5) 
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Signals  of  Type  I  and  II  are  symmetric  whereas  signals  of  Type  III  and  IV  are 
antisymmetric.  The  center  of  symmetry /antisymmetry  is  given  by  cs  =  ^f^-  +  ss. 
Notice  that  for  Type  II  and  Type  IV  signals  the  center  of  symmetry/antisymmetry 
is  in  between  two  samples. 

Given  the  above  signal  types  the  following  observations  are  in  order: 

1.  A  mirror  extended  signal  is  a  signal  of  Type  II. 

2.  The  impulse  response  of  filter  U3  f  is  a  signal  of  Type  I. 

3.  The  impulse  responses  of  filters  Ujgv  are  signals  of  Type  III. 

4.  The  circular  convolution  of  a  Type  II  signal  with  a  Type  I  filter  yields  a  Type 
II  signal. 

5.  The  circular  convolution  of  a  Type  II  signal  with  a  Type  III  filter  yields  a  Type 
IV  signal. 

As  a  consequence  of  these  results  we  have  that  the  Sj  are  always  of  Type  II  and  the 
rVj  are  always  of  Type  IV  as  long  as  the  input  signal  is  of  Type  II.  The  importance 
of  this  observation  resides  in  the  fact  that  even  though  the  mirror  extended  signal  is 
2K  samples  long,  only  K  samples  are  needed  for  computing  and  storing  the  outputs 
of  the  filter  bank. 

Now  we  are  in  a  position  to  give  a  justification  for  limiting  the  number  of 
levels  to  L  even  when  we  use  as  input  the  periodic  extension  of  a  mirror  signal.  The 
"extended"  input  signals  discussed  so  far  are  better  viewed  as  boundary  conditions 
to  constrain  the  output  signals  of  the  filter  bank  to  have  the  same  length  as  the 
input.  In  the  case  of  the  periodic  extension  it  is  clear  that  if  the  fundamental  period 
of  the  input  signal  is  shorter  than  some  filter  kernel  then  this  amounts  to  a  wrap 
around  of  the  filter  coefficients.  A  similar  result  holds  true  if  the  input  signal  is  the 
periodic  extension  of  a  mirror  signal.  If  some  filter  is  longer  than  half  the  length 
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of  the  fundamental  period  then  this  amounts  to  a  reflection  of  the  filter  coefficients 
extending  beyond  half  the  signal  length.  Since  Equation  5.5  is  based  on  the  length 
of  the  original  signal  it  is  valid  in  both  cases. 

Now  we  turn  into  a  specific  example  and  show  how  this  method  can  be  used. 
For  the  purposes  of  this  example  we  will  focus  on  a  sine-Gabor  wavelet  with  cr0  = 
4  and  lu0  =  ~- .  This  leads  to  a  narrow-band  wavelet  with  er(^)  =  0.1768  and 
relative  bandwidth  of  0.0750,  Q  =  13.33331.  In  order  to  generate  a  nearly  tight 
frame  (B/A=  1.0001)  we  need  in  this  case  at  least  8  voices  (N  =  8).  Next,  we  need 
to  sample  and  truncate  each  of  the  wavelets  ipi(t),  ...  ,  ^8(*)-  To  achieve  this  we 
extend  the  procedure  described  in  Section  5.1.  In  the  case  of  multiple  voices,  the 
mean  and  standard  deviation  of  the  sine-Gabor  wavelet  and  its  Fourier  transform  are 
given  by 

mf>„)  =  2-("-1)/iVm(</0, 

and 

=  2-("-1>/JV<r(i), 

where  tpv(t)  is  defined  by  Equation  2.35.  We  determine  the  maximum  frequency  for 
the  multivoice  sine-Gabor  wavelet  by 

where  the  value  is  such  that  the  energy  of  V;i(^)  m  the  interval  [— Umax,  Umax] 
is  larger  than  a  certain  percentage  of  its  total  energy.  Similarly,  we  truncate  the 
sine-Gabor  wavelet  ij)v{t)  to  the  time  interval 

[-dto{ipv),dto{il)v)], 

lQ  is  the  inverse  of  the  relative  bandwidth.  A  filter  is  said  to  be  narrow-band  if  Q  >  1;  broad- 
band if  Q  w  1. 
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Figure  5.6:  Magnitude  of  the  discrete  Fourier  transform  of  the  filters  gv  in  Table  5.3 
and  the  Lagrange  a  trous  filter  /  in  Table  5.4. 

where  the  value  dt  is  such  that  the  energy  of  ipv(t)  in  the  above  interval  is  larger  than 
a  certain  percentage  of  its  total  energy. 

One  may  verify  that  for  4,  =  4  more  than  99.99%  of  the  energy  of  V>i(u;)  is 
confined  to  the  interval  u)  e  [-3.07, 3.07].  Therefore  a  sampling  period  Ts  =  1  would 
be  sufficient.  Notice  that  the  bandwidth  of  wavelets  t/>2(<),  ■  •  •  ,  fo(t)  is  smaller 
than  the  bandwidth  of  ipi(t)  and  therefore  the  same  sampling  period  would  be  also 
adequate.  Next,  we  need  to  determine  the  number  of  samples  that  will  be  used  to 
define  each  gv-k  =  ^v{-k).  (Notice  that  we  have  dropped  the  complex  conjugate  from 
the  definition  of  gv.k  since  the  wavelets  in  our  study  are  real.)  As  before  we  use  an 
energy  criterion.  One  may  verify  that  for  dt  =  4  more  than  99.99%  of  the  energy  of 
ij)v{t)  is  confined  to  the  interval  t  €  [-2^-^/^11.32, 2(V-1>/Ar11.32].  Table  5.3  shows 
the  filter  coefficients  for  gv,k  for  ke[-Mv,  Mv\  where  Mv  =  [2(t,-1'/7V11.32] .  Figure 
5.6  shows  the  magnitude  of  the  discrete  Fourier  transform  of  the  filters  gv  in  Table 
5.3  and  the  Lagrange  a  trous  filter  /  in  Table  5.4.  Figure  5.7  shows  the  impulse 
response  of  the  filter  bank  in  Figure  5.5  for  the  filter  coefficients  in  Tables  5.3  and 
5.4. 
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Table  5.4:  Filter  coefficients  of  the  Lagrange  a  trous  filter  /  of  length  P=39. 
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Figure  5.7:  Response  of  the  filter  bank  in  Figure  5.5  for  the  impulse  signal  in  the  top 
left  corner  and  the  filter  coefficients  in  Tables  5.3  and  5.4. 
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5.3    Inverse  Discrete  Sine-Gabor  Wavelet  Transform 

In  this  section  we  address  technical  and  implementation  issues  regarding  the 
inverse  discrete  sine-Gabor  wavelet  transform.  A  number  of  conditions  for  the  for- 
ward transform  to  be  invertible  were  reviewed  in  Chapter  3.  We  start  by  making 
the  observation  that  for  the  forward  transform  to  be  invertible  the  discrete  Fourier 
transform  (DFT)  of  the  filters  /,  gv  cannot  be  all  zero  for  any  given  w;  for  otherwise 
the  entire  frequency  content  of  the  signal  is  not  covered  by  the  transform.  It  is  clear 
that  this  condition  is  not  satisfied  by  the  filters  in  Figure  5.6.  Therefore,  the  forward 
transform  using  these  filters  is  not  invertible.  This  might  be  a  source  of  confusion 
since  we  showed  that  the  sine-Gabor  function  is  a  wavelet  frame,  i.e.,  it  is  complete 
and  stable.  The  problem  evidenced  in  Figure  5.6  arises  from  the  disc  retization  of 
the  continuous  wavelet  transform  described  in  Chapter  3.  As  we  already  pointed  out 
in  Section  5.2,  sampling  xpv{t)  leads  to  filters  gv  which  are  of  Type  III.  It  is  easy  to 
verify  that  the  DFT  of  such  filters  is  always  zero  at  u  =  n.  Furthermore,  one  may 
verify  that  the  DFT  of  the  Lagrange  a  trous  filters  /  is  also  zero  at  w  =  jr.  One 
may  wish  to  use  an  alternative  filter  /  with  nonzero  response  at  u  =  7r.  However, 
in  order  for  the  composition  of  DWT-1  and  DWT  to  converge  rapidly,  the  response 
of  the  alternative  filter  /  at  u  =  tt  must  be  significantly  different  from  zero.  This 
defeats  the  spirit  of  the  multiresolution  approach  where  the  filters  /,  gv  split  the 
frequency  content  of  the  signal  into  low  and  high  frequencies,  respectively.  In  fact, 
the  author  is  not  aware  of  any  such  filter  /  with  nonzero  response  at  w  =  tt  leading 
to  a  reasonably  behaved  scaling  function.  We  summarize  the  above  discussion  by 
making  the  following  proposition: 

Proposition  5.3.1  The  forward  discrete  wavelet  transform  is  not  invertible  if  the 
filter  f  satisfies  f{ej7T)  =  0  and  the  filters  gv  are  defined  by  gv.k  =  ipv{-k),  k  <E 
[-M„,  Mv],  where  xpv(t)  is  any  real  odd  function,  i.e.,  tpv(t)  =  -ipv(-t),  generating  a 
wavelet  frame. 
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Proof:  The  proof  follows  from  the  above  discussion.  ■ 

Therefore,  it  follows  from  Proposition  5.3.1  that  it  is  not  possible  to  invert 
the  forward  discrete  sine-Gabor  wavelet  transform  by  directly  using  the  techniques 
reviewed  in  Chapter  3.  In  order  to  overcome  this  difficulty  we  modify  the  definition 
of  the  filters  gv  as  follows 

9v,k  =  Vv  (-(k  +  0.5)7;) ,  k  e  [-Mv,  Mv  -  1],  (5.6) 

where  Mv  =  ["2","1)gd'gWj .  This  leads  to  filters  of  Type  IV  whose  response  is 
generally  nonzero  at  u  =  7r.  Notice  that  we  have  included  the  sampling  period  as  we 
will  be  considering  samplings  other  than  Ts  =  1.  It  is  worth  pointing  out  that  the 
alternative  sampling  in  Equation  5.6  is  equivalent  to  the  sampling  shift  strategy  used 
in  [Mal92b]. 

In  order  to  understand  the  effects  of  the  sampling  shift  strategy,  one  may  view 
Equation  5.6  as  sampling  a  shifted  version  of  ipv(t),  i.e., 

9v.k  =  v*  (-kTs) ,  k  e  [-Mv,  Mv  -  1], 

where  ij>v(t)  =  ipv(t-0.5Ts).  Therefore,  it  follows  that  the  results  described  in  Chapter 
3  hold  for  i})sv{t).  Furthermore,  it  is  not  difficult  to  see  that  the  good  time-frequency 
localization  properties  of  ipv(t)  are  preserved  in  However,  two  disadvantages 

arise  from  the  alternative  sampling  in  Equation  5.6:  (1)  The  elegant  decomposition 
algorithm  described  in  Section  5.2  can  no  longer  be  used  since  the  filters  gv  are  now 
of  Type  IV,  and  (2)  The  filters  gv  are  no  longer  constant  phase.  Before  considering 
the  issues  regarding  the  inverse  transform  we  describe  how  to  successfully  evade  these 
two  problems. 
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The  first  problem  can  be  solved  by  making  the  following  two  observations: 

1.  The  circular  convolution  of  a  mirror  extended  signal  of  length  2K  (i.e.,  a  Type 
II  signal)  with  a  Type  IV  filter  of  shorter  length  yields  a  response  rk  satisfy- 
ing [Kor96] 


where  sr  is  some  integer  shift.  It  follows  from  the  above  expression  that  rST+K-X 
is  also  zero.  We  will  refer  to  the  above  response  as  a  modified  Type  III  signal. 
Notice  that  although  the  mirror  extended  signal  is  2K  samples  long,  only  K 
and  K  -  1  samples  are  needed  to  compute  and  store  the  output,  respectively. 

2.  For  a  filter  gv  of  Type  IV  one  may  show  that  U3gv  is  of  Type  III  for  j  >  0. 

From  Section  5.2  we  know  that  the  are  always  of  Type  II.  It  follows  from 
the  two  observations  above  that  the  rv.0  are  modified  Type  III  signals  and  that  the 
r„;j,  j  >  0  are  Type  IV  signals.  Therefore,  given  a  mirror  extended  signal  of  length 
2K  is  still  possible  to  have  an  elegant  algorithm  that  requires  only  K  samples  for 
computing  the  output  and  either  A'  -  1  or  K  samples  for  storing  it. 

In  order  to  understand  the  impact  of  the  second  problem,  Figure  5.8  shows 
the  step  response  of  a  3-level  decomposition  for  the  filter  coefficients  in  Table  5.3 
and  equivalent  filter  coefficients  derived  using  Equation  5.6  with  <r0  =  4,  u>o  =  -f, 
dw  =  4,  dt  =  4,  and  Ts  =  1.  It  is  clear  from  the  figure  that  by  loosing  the  constant 
phase,  the  filters  defined  by  Equation  5.6  introduce  a  shift  in  time  that  makes  tracing 
the  evolution  of  the  coefficients  in  the  time-frequency  plane  more  difficult.  However, 
this  problem  can  be  easily  evaded  by  noting  that  the  phase  of  filter  UJgv  is  given 
by  2j~1u  -  7r/2.  For  j  >  0  the  linear  part  corresponds  to  an  integer  shift  (to  the 
left)  of  2j  l.  This  can  be  compensated  for  by  shifting  the  coefficients  to  the  right 
by  the  same  amount  (see  Figure  5.8(c)).  For  j  —  0  the  shift  is  not  an  integer  and 


ke  [0,A'-1], 
k  =  2K  -  1. 
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Figure  5.8:  Step  response  of  a  3-level/8-voice  filter  bank.  The  horizontal  axis  repre- 
sents time.  The  vertical  axis  represents  the  filter  output.  From  top  to  bottom  the 
v  +  jN  row  corresponds  to  the  filter  output  rv-j,  where  N  is  the  number  of  voices. 
The  black  vertical  line  across  the  middle  of  the  image  denotes  the  location  of  the  step, 
(a)  Using  filters  from  Table  5.3.  (b)  Using  equivalent  filters  defined  by  Equation  5.6. 
(c)  Using  the  same  filters  as  in  (b)  with  shift  compensation  for  analysis  purposes  (see 
text). 
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Figure  5.9:  Magnitude  of  the  discrete  Fourier  transform  of  the  equivalent  filters  gv 
defined  by  Equation  5.6  and  the  Lagrange  a  trous  filter  /  in  Table  5.4. 

we  opt  to  leave  the  coefficients  unchanged.  Figure  5.9  shows  the  magnitude  of  the 
discrete  Fourier  transform  of  the  equivalent  filters  gv  defined  by  Equation  5.6  and 
the  Lagrange  a  trous  filter  /  in  Table  5.4. 

The  issue  of  inverting  the  transform  can  now  be  addressed  since  the  entire 
frequency  line  is  covered  by  the  filters  f,gv  (see  Figure  5.9).  Based  on  Equation 
3.52,  we  consider  the  inverse  double  integral  type  filters 

f  =  6/V2cindgv  =  cgl  (5.7) 

where  the  value  of  c  is  determined  by  finding  the  minimum  of  the  absolute  value  of 
the  difference  of  Equation  3.54  and  one  using  the  Simplex  search  method.  Notice 
that  the  factor  2(„_\)/JV  in  Equation  3.52  is  not  included  here  because  the  filters  gv  are 
defined  in  terms  of  ipv(t)  in  Equation  2.35  which  includes  an  extra  factor  u(v_1)/N 
when  compared  to  the  definition  of  tpv(t)  used  in  Section  3.4.  As  discussed  in  Section 
3.4  the  filters  in  Equation  5.7  are  only  an  approximate  inverse  for  the  discrete  wavelet 
transform.  In  order  to  evaluate  the  quality  of  the  inverse  filters  we  measure  the  error 
between  the  input,  s,  and  output,  s,  of  the  composition  of  DWT-1  and  DWT  for  a 
given  input  signal.  For  this  purpose  we  use  the  normalized  error  defined  by  [She93] 


sn  sn 
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and  evaluate  the  root  mean  square  error  given  by  [She93] 


( 


1/2 


rms 


n 


where  s  =  (||s||/||S||)s.  In  our  analysis  we  use  an  impulse  as  input  signal  and  limit 
the  number  of  levels  of  decomposition/reconstruction  to  one. 

Results  show  that  for  the  filters  discussed  above  the  mean  square  error  is  ap- 
proximately 0.4960.  The  speed  of  convergence  of  Neumann  iterations  (see  Section 
3.4.3)  is  also  relatively  slow;  the  error  after  100  iterations  with  //  =  0.5  is  approxi- 
mately 0.2947.  The  main  reason  for  this  poor  performance  is  that  the  DFT  of  the 
composition  of  DWT-1  and  DWT  for  the  filters  under  consideration  is  nearly  zero 
around  to  =  n.  The  only  way  to  improve  the  performance  is  to  sample  at  a  lower 
rate;  therefore  increasing  the  amount  of  aliasing  and  consequently  the  response  of 
the  composition  at  uo  =  n. 

Sampling  at  a  lower  rate  amounts  to  selecting  a  smaller  value  for  dw.  Filters 
given  by  Equations  5.6  and  5.7  with  a0  =  4,  uj0  =  -f,  du  =  1.5,  dt  =  4,  and 
Ts  =  1.1985  yield  a  reduced  mean  square  error  of  approximately  0.1425.  The  speed 
of  convergence  in  this  case  is  also  significantly  faster;  the  error  after  5  iterations  with 
fj.  =  0.5  is  approximately  0.0223.  The  magnitude  of  the  discrete  Fourier  transform 
of  the  above  filters  is  shown  in  Figure  5.10.  The  output  of  the  composition  and  the 
response  after  5  iterations  using  an  impulse  as  input  are  shown  in  Figure  5.11. 

Clearly,  sampling  the  wavelets  ipv{t)  at  a  lower  rate  introduces  a  certain 
amount  of  aliasing  in  the  higher  frequency  filters  gx  and  g2  (see  Figure  5.10).  Conse- 
quently, the  approximating  wavelets  ip[(t)  and  xp'2{t)  generated  by  these  filters  are  no 
longer  "close"  to  the  original  wavelets  ipi(t)  and  ip2(t)-  Figure  5.12  shows  the  discrete 
filters  gin  the  approximating  wavelets  ip'v(t),  and  the  original  wavelets  ipv{t). 

From  the  above  analysis  we  observe  that  there  is  a  trade-off  between  the  qual- 
ity of  the  inverse  filters  and  the  speed  of  convergence,  and  the  quality  of  the  higher 
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Figure  5.10:  Magnitude  of  the  discrete  Fourier  transform  of  the  equivalent  filters  gv 
defined  by  Equation  5.6  with  a0  =  4,  coQ  =  ^,  =  1.5,  dt  =  4,  and  Ts  =  1.1985, 
and  the  Lagrange  a  trous  filter  /  in  Table  5.4. 
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100  120 


(b) 

Figure  5.11:  Output  of  the  composition  of  DWT-1  and  DWT  for  the  filter  in  Figure 
5.10.  (a)  No  iterations,  (b)  After  5  iterations. 
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Figure  5.12:  Discrete  filters  gv  ('+'),  4  iterations  of  the  approximating  wavelets 
ip'v(—(t  +  0.5TS))  (red)  and  original  wavelet  ipv(  —  (t  +  0.5TS))  (blue)  with  a0  =  4, 
u0  =  ch  =  1.5,  dt  =  4,  and  Ts  =  1.1985.  From  top  to  bottom,  left  to  right,  v  =  1 
through  v  —  8. 
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frequency  approximating  wavelets.  By  improving  the  quality  of  the  inverse  filter, 
we  indirectly  change  the  shape  of  the  higher  frequency  approximating  wavelets  and 
consequently  their  time-frequency  localization  properties.  Clearly,  the  interpretation 
of  the  wavelet  coefficients  changes  for  the  affected  voices  since  the  higher  frequency 
analyzing  wavelets  are  no  longer  "close"  to  the  original  wavelets.  In  fact,  the  de- 
composition can  now  be  viewed  as  a  hybrid  wavelet  transform  in  which  the  analyzing 
wavelet  for  each  voice  is  not  necessarily  a  fractionally  dilated  version  of  a  single  mother 
wavelet.  We  remark  that  although  the  above  discussed  trade-off  and  interpretation 
are  original  to  this  work  they  also  apply  to  the  work  of  Shensa  [She93,  She96]. 


CHAPTER  6 
APPLICATIONS 

In  this  chapter  we  investigate  the  application  of  the  sine-Gabor  wavelet  trans- 
form to  image  processing.  We  describe  an  extension  of  the  transform  to  two  dimen- 
sions and  demonstrate  the  use  of  the  transform  with  both  broad-band  and  narrow- 
band sine-Gabor  wavelets  having  nearly  optimum  time-frequency  localization.  We 
also  revisit  the  use  of  the  sine-Gabor  wavelet  as  an  approximation  to  the  first  deriva- 
tive of  a  Gaussian  and  demonstrate  its  use  as  a  multivoice  multiresolution  edge  de- 
tector. 

6.1    Extension  to  Two  Dimensions 

In  order  to  extend  the  transform  to  two  dimensions  we  adopt  the  approach 
used  by  Mallat  [Mal89a]  and  construct  a  separable  two  dimensional  scaling  function 
<J>(x,  y)  =  <f>(x)(j)(y)  and  wavelets 

«(x,y)  =  rpVx(x)(f,(y),  vx  =  L...  ,N}  ,  (6.1) 

{*ly(x,y)  =  <i>(x)xl>Vv(y),vy  =  l,...,N},  (6.2) 

and 

{*lXtVv{x,  y)  =  ^(x)^(y),  vx  =  l,...,N,vy  =  l  n)  , 

where  N  is  the  number  of  voices.  For  iV  =  1  one  obtains  the  separable  extension  to 
two  dimensions  of  Mallat  [Mal89a];  that  is,  $(x,y)  =  (t>(x)(j)(y),  ^\x,  y)  =  i>(x)(f)(y), 
\f!2(x,  y)  =  4>{x)ip{y),  and  ty3(x,y)  =  ip(x)ip(y).  Notice  that  in  two  dimensions  the 
number  of  wavelets  grows  quadratically  with  respect  to  the  number  of  voices.  It  is 
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Figure  6.1:  Frequency  partition  of  a  separable  two  dimensional  scaling  function  and 
wavelets,  (a)  N  =1.  (b)  Ar  =  8. 
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Figure  6.2:  Mean  square  error  as  a  function  of  the  number  of  iterations  for  discrete 
filters  gv  with  a0  =  5/2,  w0  =  1,  dw  =  1.5,  dt  =  4,  and  Ts  =  2.2080  and  an  a  trous 
filter  /  with  Q  =  4.  From  top  to  bottom:  N  =  4,  3,  2, 1. 


Figure  6.3:  Discrete  filter  g  ('+'),  4  iterations  of  the  approximating  wavelet  ip'(-(t  + 
0.5TS))  (red)  and  original  wavelet  +  0.5TS))  (blue)  with  a0  =  5/2,  u;0  =  1, 

^  =  1.5,  dt  =  4,  and  Ts  =  2.2080. 

easy  to  verify  that  for  iV  voices  the  total  number  of  wavelets  is  N(N  +  2).  Figure  6.1 
shows  the  ideal  frequency  cover  of  a  separable  two  dimensional  scaling  function  and 
wavelets  for  TV  =  1  and  N  =  8. 

In  terms  of  its  discrete  implementation  the  two  dimensional  multivoice  wavelet 
transform  using  the  above  separable  construction  can  be  achieved  by  one  dimensional 
filtering  operations  in  the  x  and  y  directions.  Although  tedious  due  to  the  excessive 
number  of  wavelets  for  N  >  2,  the  discrete  implementation  of  the  multivoice  wavelet 
transform  is  a  trivial  extension  of  the  classical  two  dimensional  discrete  wavelet  trans- 
form [Mal89a].  To  be  concise  we  omit  the  irksome  technical  details. 

6.2    2  D  Broad-Band  and  Narrow-Band  Sine-Gabor  Wavelet 

In  this  section  we  demonstrate  the  application  of  both  broad-band  and  narrow- 
band sine-Gabor  wavelets  to  image  processing. 
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(b) 

Figure  6.4:  Four  level  decomposition  and  reconstruction  with  one  voice  (N  =  1) 
using  a  discrete  filter  g  with  a0  =  5/2,  uo0  =  1,  du  =  1.5,  dt  =  4,  and  Ts  =  2.2080  and 
an  a  trous  filter  /  with  Q  —  4.  (a)  Decomposition.  From  left  to  right:  horizontal, 
vertical,  and  diagonal  channels  were  applicable,  (b)  Reconstruction.  From  left  to 
right:  iterations  0,  1,  and  2. 
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Figure  6.5:  Mean  square  error  as  a  function  of  the  number  of  iterations  for  discrete 
filters  gv  with  a0  =  4,  uj0  =  37r/4,  dw  =  1.5,  dt  =  4,  and  Ts  =  1.1985  and  an  a  trous 
filter  /  with  Q  =  10.  From  top  to  bottom:  N  =  3  (yellow),  N  =  8  (blue),  N  =  7 
(green),  N  =  6  (red),  N  =  4  (magenta),  and  iV  =  5  (cyan). 

As  discussed  in  Chapter  4  the  sine-Gabor  wavelet  with  cr0  =  5/2  and  lo0  =  1 
yields  nearly  optimum  time-frequency  localization.  It  can  be  verified  that  Q  =  3.5355, 
i.e.,  the  wavelet  is  broad-band.  Figure  6.2  shows  the  speed  of  convergence  of  the  1- 
D  Neumann  iterations  as  a  function  of  the  number  of  voices  for  discrete  filters  gv 
with  ct0  =  5/2,  wo  =  1,  dw  =  1.5,  dt  =  4,  and  Ts  =  2.2080  and  an  a  trous  filter  / 
with  Q  =  4.  Notice  that  the  best  speed  of  convergence  is  obtained  using  only  one 
voice  (TV  =  1).  Figure  6.3  shows  the  approximating  wavelet  and  original  wavelet 
in  this  case.  Even  though  d^  is  rather  small  (approximately  93%  of  the  energy 
of  the  wavelet  is  covered  in  the  interval  [— u)max,u}max\)  the  approximating  wavelet 
is  remarkably  close  to  its  original  counterpart.  Figure  6.4  shows  a  2-D  four  level 
decomposition  and  reconstruction  using  the  wavelet  described  above.  Observe  that 
the  reconstruction  is  remarkably  good  even  without  Neumann  iterations. 

As  discussed  in  Chapter  5  the  sine-Gabor  wavelet  with  a0  =  4  and  u>o  =  3n/4 
also  yields  nearly  optimum  time-frequency  localization.  It  can  be  verified  that  Q  = 
13.3333,  i.e.,  the  wavelet  is  narrow-band.  Figure  6.5  shows  the  speed  of  convergence 
of  the  1-D  Neumann  iterations  as  a  function  of  the  number  of  voices  for  discrete 
filters  gv  with  a0  =  4,  u0  =  37r/4,  d^  =  1.5,  dt  =  4,  and  Ts  =  1.1985  and  an  a  trous 
filter  /  with  Q  =  10.  Notice  that  the  best  speed  of  convergence  is  obtained  using 
five  voices  (N  =  5).  Figure  6.6  shows  the  approximating  and  original  wavelets  in 
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Figure  6.6:  Discrete  filter  gv  ('+'),  4  iterations  of  the  approximating  wavelet  4>'v(— (t+ 
0.5TS))  (red)  and  original  wavelet  i>v{  —  {t  +  0.5TS))  (blue)  with  er0  =  4,  u>Q  =  37r/4, 
dw  =  1.5,  dt  =  4,  and  Ts  =  1.1985.  From  top  to  bottom,  left  to  right,  v  =  1  through 

v  =  5. 
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Figure  6.7:  Two  level  decomposition  and  reconstruction  with  5  voices  (N  =  5)  using 
discrete  filters  gv  with  a0  =  4,  u0  =  37r/4,  =  1.5,  dt  =  4,  and  Ts  =  1.1985  and  an 
a  trous  filter  /  with  Q  =  10.  (a)  Input  image  (2-D  chirp),  (b)  Second  level  channels. 
The  horizontal  and  vertical  labels  (dcx,  vx  =  5,4,  and  dcy,  vy  =  5,4)  indicate  the 
channel  shown. 
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Figure  6.7:  Continued:  (c)  Second  level  channels.  The  horizontal  and  vertical  labels 
(vx  =  3,  2, 1,  and  dcy,  vy  =  5,  4)  indicate  the  channel  shown. 
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Figure  6.7:  Continued:  (d)  Second  level  channels.  The  horizontal  and  vertical  labels 
(dcx,  vx  =  5,  4,  and  vy  =  3,  2, 1)  indicate  the  channel  shown. 
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Figure  6.7:  Continued:  (e)  Second  level  channels.  The  horizontal  and  vertical  labels 
(vx  =  3,  2, 1,  and  vy  =  3, 2, 1)  indicate  the  channel  shown,  (f)  Reconstruction.  From 
left  to  right:  iterations  0,  2,  and  5. 
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Figure  6.8:  Mean  square  error  as  a  function  of  the  number  of  iterations  for  discrete 
filters  gv  with  a0  =  1.0657,  lj0  =  0.0299,  rfw  =  1,  dt  =  4,  and  Ts  =  2.0862  and  an  a 
trous  filter  /  with  Q  =  2.  From  top  to  bottom:  N  =  3, 2, 1. 

this  case.  Observe  that  for  v  >  3  the  approximating  wavelets  are  remarkably  close 
to  their  original  counterparts.  Figure  6.7  shows  a  2-D  two  level  decomposition  and 
reconstruction  using  the  wavelet  described  above. 

Based  on  the  above  examples  and  other  experimental  data  we  observed  that: 

1.  Broad-band  sine-Gabor  wavelets  are  more  suitable  for  applications  were  good 
time  resolution  is  required.  Usually  only  one  voice  is  required;  although  multiple 
voices  can  be  used  to  approximate  a  continuous  transform  at  the  expense  of 
slower  convergence.  Figure  6.4  suggests  broad-band  sine-Gabor  wavelets  may 
be  used  for  boundary  detection. 

2.  Narrow-band  sine-Gabor  wavelets  are  more  suitable  for  applications  were  good 
frequency  resolution  is  required.  Usually  multiple  voices  are  required  to  obtain 
good  convergence.  Figure  6.7  suggests  narrow-band  sine-Gabor  wavelets  may 
be  used  for  texture  segmentation/classification. 

3.  By  modifying  the  parameters  a0  and  u)0  it  is  possible  to  change  the  relative 
bandwidth  (or  Q)  of  the  sine-Gabor  wavelet  transform  while  preserving  nearly 
optimum  time- frequency  localization. 
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Figure  6.9:  Discrete  filter  gv  ('+'),  4  iterations  of  the  approximating  wavelet  ip'v{—(t+ 
0.5TS))  (red)  and  original  wavelet  ijjv(—(t  +  0.5TS))  (blue)  with  cr0  =  1-0657,  u0  = 
0.0299,  dw  =  1,  dt  =  4,  and  Ts  =  2.0862.  From  top  to  bottom  v  =  1  through  v  =  3. 

6.3    2  D  Sine  Gabor  Edge  Detector 

In  this  section  we  investigate  the  use  of  the  sine-Gabor  wavelet  as  an  edge 
detector.  Canny  [Can86]  introduced  an  optimal  detector  for  1-D  step  edges  and  sug- 
gested the  first  derivative  of  a  Gaussian  as  an  efficient  approximation.  Image  edges 
of  a  particular  orientation  can  be  detected  by  convolving  the  image  with  a  mask 
generated  by  convolving  a  linear  edge  detection  function  aligned  normal  to  the  edge 
direction  with  a  projection  function  parallel  to  the  edge  direction1  [Can86].  Signifi- 
cant computational  savings  are  possible  if  the  projection  function  is  a  Gaussian  with 
the  same  a0  as  the  first  derivative  of  the  Gaussian  used  as  the  detection  function 
[Can86].  In  this  case  edges  of  a  particular  orientation  can  be  detected  by  convolv- 
ing the  image  with  a  symmetric  two  dimensional  Gaussian  and  then  differentiating 
normal  to  the  edge  direction.  In  fact  this  procedure  need  not  be  carried  out  in  every 
direction  because  the  slope  of  a  smooth  surface  in  any  direction  can  be  determined 
exactly  from  its  slope  in  two  directions  [Can86].  The  operation  can  be  summarized 
as  follows  [Can86,  Mal92b]: 

1.  An  image  s(x,y)  is  smoothed  via  convolution  with 

e(x,y)  =  KGexp(-^f\.  (6.3) 

lrThe  term  edge  direction  refers  to  the  direction  of  the  tangent  to  the  contour  that  the  edge 
defines  in  two  dimensions. 
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Figure  6.10:  Three  level  decomposition  and  reconstruction  with  3  voices  (N  =  3) 
using  discrete  filters  gv  with  a0  =  1.0657,  lo0  =  0.0299,  du  =  1,  dt  =  4,  and  Ts  = 
2.0862  and  an  a  trous  filter  /  with  Q  =  2.  The  horizontal  and  vertical  labels  (dcx, 
vx  =  3,2,1,  and  dcy,  vy  =  3.2,1)  indicate  the  channel  shown,  (a)  Input  image, 
(b)  Level  1. 
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Figure  6.10:  Continued:  (c)  Level  2. 


Figure  6.10:  Continued:  (d)  Level  3.  (e)  Reconstruction.  From  left  to  right:  iterations 
0,  2,  and  5. 
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Figure  6.12:  Gradient  magnitude,  gradient  angle,  and  edges  of  a  three  level  de- 
composition with  three  voices  (N  =  3)  using  discrete  filters  gv  with  a0  =  1.0657, 
lo0  =  0.0299,  dw  =  1,  dt  =  4,  and  Ts  =  2.0862  and  an  a  trous  filter  /  with  Q  =  2. 
From  left  to  right:  magnitude,  angle,  and  edges.  From  top  to  bottom:  v  =  1,2,3. 
The  angle  is  shown  only  when  the  magnitude  is  significantly  different  from  zero.  An- 
gles within  135  ±  22.5  degrees,  90  ±  22.5  degrees,  45  ±  22.5  degrees,  and  0  ±  22.5 
degrees  are  shown  in  yellow,  red,  cyan,  and  blue,  respectively,  (a)  Level  1. 
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Figure  6.12:  Continued:  (c)  Level  3. 
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2.  The  gradient  vector  is  then  computed 


V(s*6){x,y) 


The  direction  of  the  gradient  vector  at  a  point  (x0,  y0)  indicates  the  direction  in 
the  image  plane  along  which  the  directional  derivative  of  f{x,  y)  has  the  largest 
absolute  value. 

3.  Edges  are  defined  as  points  (x0,  yQ)  where  the  modulus  of  the  gradient  vector  is 
maximum  in  the  direction  of  the  gradient  vector;  i.e.,  edge  points  are  inflection 
points  of  the  surface  (s  *  9)(x,  y). 

The  procedure  is  usually  carried  out  at  different  scales  to  find  the  best  trade-off 
between  detection  and  localization  [Can86]. 

Mallat  [Mal92b]  related  the  above  algorithm  to  the  formulation  of  a  2-D 
wavelet  transform  by  defining  two  wavelet  functions  ^(x,y)  and  ^2(x,y)  such  that 


The  wavelet  decomposition  using  these  wavelets  is  equivalent  to  a  multiscale  Canny 
edge  detector  [Mal92b].  Due  to  the  separability  of  9(x,  y)  in  Equation  6.3  the  wavelets 
^(x,  y)  and  *2(x,  y)  may  be  written  as  in  Equations  6.1  and  6.2  where  0  and  ip  are 
a  Gaussian  and  its  first  derivative,  respectively. 

As  discussed  in  Chapter  4  the  sine-Gabor  wavelet  with  cr0  =  1.0657  and  ujq  = 
0.0299  is  for  practical  purposes  equal  to  the  first  derivative  of  a  Gaussian  (Canny's 
approximation  to  the  optimal  step  edge  detector  [Can86]).  Figure  6.8  shows  the  speed 
of  convergence  of  the  1-D  Neumann  iterations  as  a  function  of  the  number  of  voices 
for  discrete  filters  gv  with  a0  =  1.0657,  u0  =  0.0299,  4  =  1,  dt  =  4,  and  Ts  =  2.0862 
and  an  a  trous  filter  /  with  Q  =  2.  Notice  that  the  best  speed  of  convergence  is 
obtained  using  only  one  voice  (N  =  1).  However,  in  order  to  demonstrate  the  use 
of  voices  we  select  iV  =  3  at  the  expense  of  having  slower  convergence.  Figure  6.9 


y\x,y) 


d6(x,y) 


dx 


and  *2(x,  y) 


dO(x,y) 
dy 
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shows  the  approximating  and  original  wavelets  in  this  case.  Figure  6.10  shows  a  2-D 
three  level  decomposition  and  reconstruction  using  these  wavelets.  Although  0  is 
not  a  Gaussian  we  assume  the  approximation  as  valid  (see  Figure  6.11).  Figure  6.12 
shows  the  gradient  magnitude,  gradient  angle,  and  edges  using  the  above  filters  for 
the  2-D  three  level  decomposition  shown  in  Figure  6.10.  The  resulting  edge  set  may 
be  viewed  as  the  output  of  a  multivoice  multiresolution  Canny  edge  detector.  The 
addition  of  voices  may  prove  useful  when  trancing  the  evolution  of  edges  through 
scale-space  is  required  (e.g.,  singularity  detection  and  processing  [Mal92a]). 


CHAPTER  7 
CONCLUSIONS 

In  this  work  we  introduced  a  new  nonorthogonal  wavelet:  the  sine-Gabor 
function.  We  showed  that  this  function  not  only  satisfies  the  wavelet  admissibility 
condition  but  achieves  nearly  optimum  time-frequency  localization  and  generates 
a  frame  of  wavelets.  Furthermore,  we  showed  that  there  is  a  trade-off  between  the 
tightness  of  the  frame  and  the  time-frequency  localization  properties  of  the  sine-Gabor 
wavelet.  We  demonstrated  how  this  trade-off  can  be  overcome  by  introducing  voices. 
We  also  showed  that  under  certain  conditions  the  sine-Gabor  wavelet  is  nearly  equal 
to  the  first  derivative  of  a  Gaussian  and  generates  a  "snug"  frame.  We  addressed  the 
implementation  issues  regarding  the  computation  of  the  discrete  wavelet  transform. 
We  showed  that  there  exist  a  class  of  wavelets  (the  sine-Gabor  wavelet  being  one 
of  them)  for  which  the  discrete  wavelet  transform  is  not  invertible.  We  described 
how  to  evade  this  problem  by  modifying  the  strategy  for  selecting  the  discrete  filters 
implementing  the  transform.  Furthermore,  we  showed  how  in  the  case  of  the  sine- 
Gabor  wavelet,  desirable  properties  such  as  constant  phase  and  nearly  optimum  time- 
frequency  localization  are  preserved  by  the  alternative  filters.  We  concluded  our 
investigation  by  extending  the  transform  to  two  dimensions  and  demonstrating  some 
of  its  uses  in  image  processing  including  a  multivoice  multiresolution  edge  detector. 
However,  a  number  of  topics  remain  to  be  investigated.  Understanding  the  link 
between  the  role  of  voices  in  the  continuous  and  discrete  implementations  would 
be  desirable  for  filter  design  purposes.  We  have  seen  that  although  a  wavelet  may 
generate  a  "snug"  frame,  its  discrete  implementation  is  not  necessarily  invertible  or 
converges  rapidly.  In  particular,  the  relation  between  the  speed  of  convergence  of 


108 


109 


the  discrete  implementation  and  the  tightness  of  the  wavelet  frame  needs  further 
study.  The  author  speculates  that  the  role  of  voices  in  the  discrete  case  is  somewhat 
different  than  in  the  continuous  case.  Although  in  the  continuous  case  the  addition  of 
voices  leads  to  tighter  frames,  in  the  discrete  case  this  is  not  necessarily  the  case.  We 
conjecture  that  using  a  transform  that  suppresses  a  few  low  frequency  voices  would 
yield  faster  convergence  than  a  transform  that  uses  fewer  voices  to  avoid  the  overlap 
between  the  lowpass  filter  /  and  the  low  frequency  voices. 

Other  areas  of  research  worth  pursuing  include  the  development  of  an  efficient 
algorithm  for  evaluating  Equation  3.31  given  a  signal  s(t)  or  its  sampled  version 
according  to  Shannon's  sampling  theorem,  and  studying  the  quality  of  the  approx- 
imation of  the  discrete  wavelet  transform  to  the  exact  sampled  wavelet  transform. 
Finally,  the  extension  of  the  proposed  methodology  for  2  D  nonseparable  nonorthog- 
onal  wavelets  would  certainly  be  fruitful  for  image  processing  and  computer  vision 
applications  that  use  ad-hoc  procedures  for  computing  and  inverting  2  D  nonsepa- 
rable nonorthogonal  wavelet  transforms. 
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